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In  perturbation  theory  studies  of  the  large-scale  structure  of  the  imiverse,  it 

is  common  procedure  to  ignore  the  effects  of  an  initial  rotational  velocity  field. 

This  has  been  justified  by  the  fact  that  in  the  linear  regime,  rotational  velocities 

decay  with  the  expansion  of  the  universe,  a  simple  result  of  conservation  of  angular 

momentum.  In  this  dissertation,  we  explore  the  consequences  of  retaining  this  mode, 

based  on  the  assumption  that  gravitational  collapse  should  lead  to  the  return  of 

appreciable  rotational  velocities.    We  extend  the  classical  perturbation  theory  to 

include  rotational  modes  up  to  the  quasi-hnear  third  order  and  find  that  these  modes 

contribute  in  a  non-trivial  manner.  We  then  include  these  modes  in  the  non-Hnear 

regime,  incorporating  them  into  the  halo  model.  We  study  a  specific  appUcation  of 

non-Unear  collapse,  the  calculation  of  the  kinetic  Sunyaev-Zeldovich  effect,  with  the 

inclusion  of  rotational  velocity  modes.  We  find  that  these  modes  do  not  contribute 

significantly  at  any  angular  scale  available  to  modern  observation. 


IX 


CHAPTER  1 
INTRODUCTION 

1.1     Overview 

Cosmology,  the  study  of  the  nature  and  evolution  of  the  universe  as  a  whole, 
has  emerged  in  the  last  one  hundred  years  as  a  mature  science.  Though  cosmol- 
ogy encompasses  a  great  many  fields  of  study,  we  will  primarily  be  interested  in 
what  is  called  large-scale  structure,  the  study  of  the  general  distribution  of  matter 
and  energy  in  the  universe.  To  understand  this  distribution  today,  we  must  first 
understand  what  processes  underlie  its  evolution. 

The  first  part  of  the  twentieth  century  was  marked  by  a  great  number  of  ad- 
vancements in  the  physical  sciences,  not  the  least  of  which  was  Albert  Einstein's 
development  of  the  theory  of  relativity.  This  is  a  theory  in  which  time  and  space 
are  no  longer  absolute,  and  the  physics  of  Newton  and  his  followers  no  longer  ap- 
phes.  One  of  the  strangest  predictions  of  the  theory  was  that  the  universe  may  not 
be  static,  that  is,  that  the  universe  as  a  whole  may  be  expanding  or  contracting. 
If  it  were  expanding,  then  at  some  point  in  the  past,  everything  must  have  been 
much  closer  together,  a  tiny  speck  that  expanded  into  everjdhing  in  the  universe 
today.  This  state  of  affairs  became  known  as  the  big  bang.  Einstein  himself  held 
this  possibility  in  distaste,  and  when  it  was  shown  that  a  static  universe  was,  in 
fact,  necessarily  unstable,  he  added  another  term  into  his  equations,  the  cosmolog- 
ical  constant,  to  allow  for  a  stable  universe.  Not  long  after,  the  astronomer  Edwin 
Hubble  made  the  discovery  that  galaxies  far  away  from  our  own  are  receding  from  us, 
in  fact  the  farther  away  they  are,  the  faster  they  are  receding.  The  only  conclusion 
to  be  drawn  was  that  the  universe  is  expanding.  So  the  distribution  of  matter  and 
energy  in  the  universe  is  first  influenced  by  the  general  expansion  of  the  universe. 


The  other  process  which  is  beUeved  to  play  a  major  role  in  shaping  today's 
universe  is  gravity.  It  is  believed  that  tiny  irregularities  early  in  the  universe  grew, 
as  the  universe  expanded,  by  gravitational  attraction.  A  volume  with  a  slightly 
higher  density  than  the  surrounding  region  would  grow  by  gravitational  accretion 
until  it  had  gathered  enough  to  form  the  stars  and  galaxies  we  see  today.  In  actuality, 
however,  what  we  see  is  but  a  tiny  fraction  of  all  that  is  out  there.  The  universe  is 
dominated  by  mysterious  "dark  matter"  and  even  more  mysterious  "dark  energy." 
Dark  matter  is  matter  that  we  cannot  see,  but  it  is  influenced  by  gravity  in  the  same 
way  as  normal  matter.  One  of  the  primary  goals  of  cosmology  today  is  to  determine 
the  nature  of  dark  matter.  Dark  energy  is  even  less  well  understood,  but  its  presence 
is  implied  by  a  number  of  current  intriguing  observational  and  theoretical  findings. 
It  is  believed  responsible  for  the  fact  that  the  rate  of  expansion  of  the  universe  seems 
to  be  increasing.  Curiously,  it  can  be  seen  as  a  realization  of  the  very  cosmological 
constant  that  Einstein  called  his  greatest  blunder.  ,;    \ 

With  these  two  ingredients  we  are  in  a  position  to  describe  the  evolution  of 
both  fluctuations  in  the  density  field  of  the  universe  and  fluctuations  in  the  peculiar 
velocity,  the  quantity  of  central  interest  in  this  dissertation.  The  peculiar  velocity 
is  simply  the  velocity  of  some  cosmological  object  with  the  velocity  due  to  the 
expansion  of  the  universe  factored  out. 

One  way  to  get  a  handle  on  the  nature  of  the  distribution  of  matter  and  energy 
today,  as  well  as  far  into  the  past,  is  to  look  at  fluctuations  of  the  temperature  of  the 
cosmic  microwave  background  (CMB).  The  CMB  is  a  steady  stream  of  microwave 
light  that  reaches  us  after  a  journey  of  more  than  13  billion  years.  The  light  comes 
from  a  period  of  time  called  recombination,  just  a  few  hundred  thousand  years  after 
the  big  bang,  so  it  is  able  to  indirectly  give  us  information  about  the  state  of  the 
universe  at  that  time.  It  also  contains  information  about  the  imiverse  of  the  more 
recent  past  and  present  because  some  of  it  has  interacted  with  galaxies  and  clusters 


of  galaxies  on  its  way  to  us.  The  peculiar  velocities  of  these  foreground  objects 
cause  distinct  fluctuations  in  the  temperature  of  the  CMB,  a  process  known  as  the 
kinetic  Sunyaev-Zeldovich  (kSZ)  effect.  The  particular  goal  of  this  dissertation  is  to 
determine  the  role  rotational  peculiar  velocities  play  in  this  effect.  It  is  hoped  that 
near-future  astronomical  missions  can  provide  observations  of  these  fluctuations  to 
verify  the  predictions  made  by  the  powerful  theoretical  framework. 

1.2     Motivation  and  Outline 

The  initial  conditions  for  formation  of  large-scale  structure  are  imprinted  at 
the  time  of  the  decoupling  of  matter  from  radiation.  The  peculiar  velocity  field  of 
matter  at  that  time  can  be  decomposed  into  two  modes.  The  first,  called  the  lon- 
gitudinal mode,  is  a  purely  radial  velocity.  The  other,  called  the  rotational  mode, 
is,  as  the  name  suggests,  a  pure  rotation.  In  most  cosmological  apphcations  it  has 
been  common  procedure  to  disregard  the  rotational  mode  of  the  pecuhar  velocity 
because  it  is  seen  to  decay  as  the  universe  expands,  a  simple  consequence  of  the  con- 
servation of  angular  momentum.  While  this  is  indeed  true,  conservation  of  angular 
momentum  demands  the  return  of  appreciable  rotational  velocities  upon  gravita- 
tional collapse  of  cosmological  objects.  This  means  that  for  apphcations  where 
nonlinear  (large  density  fluctuations)  scales  are  important,  the  rotational  velocities 
may  make  important  contributions.  One  application  where  this  may  be  important 
is  in  the  aforementioned  kSZ  effect,  where  photons  from  the  CMB  are  scattered  by 
moving,  collapsed,  foreground  objects.  It  is  hoped  that  the  next  generation  of  CMB 
anisotropy  observations  will  be  able  to  see  this  effect. 

The  basic  model  for  understanding  the  evolution  of  densities  and  velocities,  the 
gravitational  instability  picture,  will  be  discussed  in  chapter  2.  Chapter  3  will  ex- 
plore the  perturbation  theory  approach  with  the  rotational  velocity  modes  included. 
In  chapter  4  we  will  introduce  the  physics  of  the  CMB  and  its  anisotropies.  Chap- 
ter 5  will  discuss  the  ingredients  that  go  into  a  calculation  of  the  expected  angular 


power  spectrum  of  CMB  temperature  fluctuations  and  then  show  the  results  for 
the  quasi-hnear  and  nonhnear  regimes,  both  without  inclusion  of  rotational  velocity 
effects.  In  chapter  6  we  will  build  a  model  which  includes  rotational  velocity  effects 
and  calculate  the  expected  angular  spectrum  with  this  inclusion.  Chapter  7  will 
contain  concluding  remarks. 


CHAPTER  2 
FOUNDATIONS 

2.1     Gravitational  Instability 

One  of  the  great  cosmological  puzzles  of  the  last  half  century  has  been  ex- 
plaining the  nature  and  evolution  of  inhomogeneities  in  today's  universe.  Einstein's 
cosmological  principle  makes  the  claim  that  on  the  largest  scales  the  universe  is 
homogeneous  and  isotropic,  and  this  assumption  is  strongly  supported  by  numerous 
observations.^  On  smaller  scales,  however,  the  universe  is  clearly  quite  inhomoge- 
neous,  as  matter  is  clustered  into  stars,  galaxies,  clusters  of  galaxies,  and  beyond, 
with  structure  extending  to  the  great  walls  and  voids  seen  on  scales  of  a  few  percent 
of  the  observable  universe. 

The  leading  picture  for  understanding  the  growth  of  these  inhomogeneities  from 
the  early  universe  to  today  is  known  quite  broadly  as  the  gravitational  instability 
model.  The  basic  assumption  of  this  model  is  that  gravitation  alone  is  responsible  for 
the  growth  of  inhomogeneities  from  their  seeds  in  the  early  universe  to  their  present 
state.  The  exact  origin  of  these  seeds  is  still  something  of  an  open  question,  but 
it  is  fairly  well  accepted  today  that  they  are  inflationary  amplifications  of  quantum 
zero-point  fluctuations.^  Quite  generally,  perturbations  can  be  of  two  orthogonal 
modes,  adiabatic  or  isocurvature,  or  a  linear  combination  of  the  two.  Adiabatic  per- 
turbations are  those  where  volume  elements  are  uniformly  contracted  or  expanded, 
causing  equal  changes  to  the  number  densities  of  matter  and  radiation.  This  obvi- 
ously causes  perturbations  in  the  energy  densities,  and,  in  fact,  affects  the  energy 
density  of  matter  differently  than  that  of  radiation.  The  orthogonal  mode,  consist- 
ing of  isocurvature  perturbations,  describes  entropy  density  fluctuations  where  the 


energy  densities  remain  smooth.  Number  density  fluctuations  in  one  species  are  bal- 
anced by  those  of  the  others.  At  times  well  earlier  than  the  era  of  matter-radiation 
equivalence,  both  modes  are  distinct  and  take  on  important  roles.  At  and  after  the 
matter-radiation  equivalence  era,  we  are  principally  interested  in  sub-horizon  sized 
fluctuations,  and  sub-horizon  sized  isocurvature  perturbations  become  adiabatic  as 
pressure  differences  lead  to  actual  density  fluctuations.  Thus  beginning  at  the  era  of 
equivalence  we  can  effectively  treat  all  perturbations  as  adiabatic  and  use  Newtonian 
fluid  dynamics  to  follow  their  evolution. 

The  use  of  Newtonian  mechanics  in  place  of  general  relativity  in  describing 
the  evolution  of  density  fluctuations  is  justified  so  long  as  the  characteristic  size 
of  the  perturbations  is  small  compared  to  the  horizon  size  and  the  magnitude  of 
the  density  perturbations  is  not  large,  thus  ensuring  nonrelativistic  gravitational 
potentials.  To  find  the  equations  that  govern  such  evolution  we  begin  by  writing 
down  the  basic  equations  for  an  ideal,  pressureless  fluid.  The  first  is  the  continuity 
equation,  expressing  conservation  of  mass. 

^]    +V..(p«)  =  0,  (2.1) 


or 

d 


at 


p  =  -pVr  •  u.  (2.2) 


Here  p  is  the  density,  r  is  the  proper  distance,  and  u  its  time  derivative.  We  can 
also  write  down  the  Euler  equation,  expressing  conservation  of  momentum, 

Here  p  refers  to  pressure,  which  we  will  set  to  zero  for  the  remainder  of  this  discus- 
sion. When  density  contrasts  become  large,  as  in  the  collapse  of  a  cluster  of  material, 
pressure  again  becomes  important,  but  at  that  point  the  perturbative  approach  we 
axe  now  exploring  fails.  The  $  in  the  above  equation  is  the  gravitational  potential. 


The  final  necessary  fluid  flow  equation  is  the  Poisson  equation 

V^$  =  iirGp,  (2.4) 

where  G  is  Newton's  gravitational  constant.  Now  we  wish  to  introduce  a  change  of 
variables  to  "co- moving"  coordinates,  x  =  r/a{t).  a{t)  is  just  the  cosmological  scale 
factor,  as  it  appears  in  the  most  general  possible  line  element  for  a  homogeneous, 
isotropic,  expanding  universe,  the  Robertson-Walker  Une  element^ 

dx'' 


ds^  =  -c^dr  +  a\t) 


+  x^dQ' 


(2.5) 


(1  -  kx^) 

The  change  of  variables  to  co-moving  coordinates  is  essentially  a  "subtracting  out" 
of  the  effects  of  the  expansion  of  the  universe.  With  this  definition  of  x  we  find 

u  =  r    =    dx  +  xa 

=    Hr  +  v,  (2.6) 

where  H  =  a/a.  H  is  known  as  the  Hubble  constant,  which  is  a  function  of  time. 
It  is  normally  expressed  as  H  =  100/ikm  •  s~^Mpc~^,  with  the  present  value  of  the 
dimensionless  constant  being  h  =  0.71  ±0.04.^  The  velocity  in  these  coordinates,  v, 
is  called  the  peculiar  velocity,  and  is  just  the  velocity  of  the  material  in  excess  of  its 
expansion  velocity.  It  is  the  "perturbation"  to  the  background  smooth  expansion. 

Before  continuing  our  look  at  the  fluid  flow  equations,  we  should  pause  to  take 
a  closer  look  at  the  behavior  of  a{t)  and  its  time  derivatives.  To  begin  we  express  the 
pressureless  Euler  equation  in  terms  of  the  new  co-moving  coordinates,  noting  that 
in  the  unperturbed  case,  dx/dt  =  0,  and  that  the  solution  for  $o  in  the  unperturbed 
Poisson  equation  is  $o  =  (2/3)7rGpo^^-  This  gives  us 

^(acc)  =  ---TrGpooV,  (2.7) 


which  leaves,  (as  V(a:^)  =  2j;), 


4 
a  =  --TrGpoa.  (2.8) 


Here  the  V  refers  to  differentiation  with  respect  to  x.  We  now  have  the  behavior 
of  o,  which  we  wiU  need  later  in  this  discussion.  We  gain  further  information  by 
integrating  the  above  equation  once  over  time.  This  gives 


a\^      8    ^         C 


,    =-nGpo  +  ^.  (2.9) 

The  constant  of  integration,  C,  is  seen  to  be  proportional  to  the  k  in  the  Robertson- 
Walker  metric,  which  is  a  measure  of  the  overall  curvature  of  the  universe.  The  most 
recent  findings'*  suggest  that  the  imi verse  is  flat,  so  we  set  C  =  A;  =  0  and  ignore 
the  contribution  of  curvature  to  the  overall  energy  density  of  the  imiverse  for  the 
remainder  of  the  discussion. 

Now,  back  to  the  discussion  of  the  fluid  flow  equations.  Since  we  are  interested 
in  examining  the  evolution  of  the  density  fluctuations,  we  write  the  density  as  p  = 
Po[l  +  (5]-  Our  goal  then  is  to  determine  the  evolution  of  5  and  v.  po  is  the  mean 
density,  which  satisfies  po  oc  a{t)'^  and  therefore  has  a  time  evolution  po  =  SpoH. 
Finally,  we  note  the  identity  that  for  any  function  /  under  this  change  of  coordinates 

Now  we  are  in  a  position  to  rewrite  the  fluid  flow  equations  in  terms  of  the 

co-moving  variables  and  the  perturbations,  5  and  v.  The  continuity  equation  now 

looks  like 

(^-~x-  v)  po[l  +  S]  +  ^V-  [(1  +  S)iax  +  v)]  =  0,  (2.11) 

\at      a         J  a 


which  simphfies  to 


S+-V-[il  +  5)v]  =  0.  (2.12) 

a 


Poisson's  equation  is  now 

1 


a 


V2$  =  47rGpo(l  +  <5).  (2.13) 


^^2 

2  *' 


To  remove  the  unperturbed  part  we  write 

$  =  0  + -TrGpoaV,  (2.14) 

o 

where  the  second  term  on  the  right  hand  side  is  the  solution  to  the  unperturbed 
Poisson's  equation,  as  above.  Using  this  for  $  we  have 

VV  =  AnGpoa^S.  (2.15) 

At  this  point  we  should  note  that  in  the  above  discussion  we  have  assumed  the 
density  to  consist  only  of  mass  density.  This  is  certainly  not  at  all  the  case,  as  fully 
seventy  percent  of  the  imiverse's  energy  density  is  contained  in  a  mysterious  "dark 
energy"  or  cosmological  constant  that  provides  a  negative  pressure  and  is  thought  to 
be  fueling  the  acceleration  of  the  expansion  of  the  universe.  As  the  name  suggests, 
however,  the  contribution  to  the  energy  density  from  the  cosmological  constant  is 
smooth  and  does  not  influence  the  evolution  of  density  fluctuations;  thus  the  terms 
involving  the  cosmological  constant  that  belong  in  the  full  fluid  flow  equations  cancel 
when  we  reduce  to  the  perturbation  evolution  equations. 

Finally,  we  look  at  the  Euler  equation  under  the  change  of  coordinates.  With 
the  pressure  term  and  the  cosmological  constant  term  left  out  we  have 

I  —  -  -X  ■  v]  {dx  +  v)  +  -{dx  +  v)  ■  V(dx  +  v)  =  — V0 nGpoO^x.  (2.16) 

\<7t      a  J  a  a  a3 

Using  the  result  of  equation  [2.8]  for  the  second  time  derivative  oia{t),  this  simplifies 

to 

5v      a         1  ,      „,  1  „  ,  ,        , 

^r-  +  -v  +  -{v-V)v  =  — V0.  (2.17 

at       a         a  a 
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This  equation,  along  with  eqns.[2.12]  and  [2.15]  describe  the  evolution  of  6  and  v. 
By  combining  equations  [2.12]  and  [2.17]  we  can  write  one  second  order  differential 
equation  for  the  evolution  of  the  density  contrast. 

W  ^  '^Idi  '  ^''^'""'^  "  47rGpo<5'  +  -^V(f>  ■  VS  +  ^ViV,[(l  +  S)v'v^].      (2.18) 

This  equation  is  written  so  that  the  left  hand  side  has  only  terms  of  first  order 
in  the  perturbed  quantities,  while  the  right  hand  side  has  the  higher  order  terms. 
Assuming  the  perturbations  are  small  (i.e.  6  <^  1)  we  can,  to  good  approximation, 
set  the  right  hand  side  equal  to  zero.  When  S  approaches  unity,  the  above  equation 
can  be  used  to  generate  contributions  to  second  order  or  higher.  Analyses  of  these 
terms  fall  under  the  heading  of  "weakly  non-linear"  calculations  and  have  been  often 
explored.^^  At  <5  >  1  perturbation  theory  breaks  down  and  different  methodologies 
must  be  pursued. 

In  order  to  solve  for  the  first  order  behavior  of  S  we  need  to  go  back  and  look 
at  the  time  dependence  of  a{t).  Equation  [2.9]  provides  a  differential  equation  for 
a{t)  which  we  can  solve,  noting  again  that  po  oc  a~^.  The  solution  is  a{t)  oc  ^^/^,  so 
that  a/a  =  H  =  2/3t.  Our  differential  equation  for  S  is  now,  to  first  order, 

d^S       4d6        2  ^      ^  ,        ^ 

The  two  solutions  are  6  oc  t^^^,  the  growing  mode,  and  S  oct~^,  the  decaying  mode. 
We  will  take  a  closer  look  at  these  solutions  in  the  next  chapter,  as  well  as  solutions 
to  higher  order  in  perturbation  theory. 

With  these  solutions  in  hand  we  can  now  begin  to  look  at  the  behavior  of  the 
quantity  of  central  interest  in  this  dissertation,  the  peculiar  velocity.  We  begin  by 
noting  that  by  observation  we  can  write  the  general  solution  to  the  linearized  version 
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of  equation  [2.12]  (ignoring  the  V  •  {5v)  term)  as 

v{x,  t)  =  -aVV-^S{x,  t)  +  ^^,  (2.20) 

where  F  is  an  arbitrary  function  which  satisfies  V  ■  F  =  0.  We  will  pursue  the 
higher  order  terms  in  the  next  chapter.  Here  the  V~^  operator  simply  refers  to 
the  solution  of  the  Poisson  equation,  V~^(5  —  (p/AirGpoO^.  It  is  the  second  term  in 
this  expression  for  the  velocity  which  will  be  of  central  importance  for  the  majority 
of  this  dissertation.  We  can  see  that  it  decays  as  o~\  and  thus  is  typically  (and 
often  justifiably,  for  the  large  scales  we  are  dealing  with  in  the  perturbation  theory 
regime)  ignored.  The  decay  of  this  term  is  expected  as  a  result  of  simple  conserva- 
tion of  angular  momentum  in  an  expanding  universe.  Length  grows  by  a  factor  of 
a{t),  so  the  product  of  length  and  rotational  velocity  should  remain  constant.  See 
Appendix  A  for  a  more  rigorous  derivation  of  conservation  of  angular  momentum  in 
an  expanding  universe.  Of  course,  when  objects  collapse,  on  small  scales,  the  same 
conservation  of  momentum  demands  the  return  of  appreciable  rotational  velocities. 
We  can  write  this  solution  more  explicitly  by  writing  S{x,t)  oc  D{t)  where 
D  could  be  either  the  growing  or  decaying  time  dependence,  and  defining  /  = 
{a/D){dD/da).  If  we  only  retain  the  growing  mode  for  D,  and  we  can  ignore  the 
relativistic  background,  as  we  already  have,  then  /  can  be  approximated  by  a  simple 
function  of  Q,  /(r2)  =  ^^■^}  Additionally,  the  Poisson  equation  is  familiar  from, 
among  many  physical  applications,  electromagnetism,  and  its  solution  is 


J  \x' 


ct>  =  -poCa^  /  rfV^^^.  (2.21) 


So,  ignoring  the  divergence-free  part  of  the  velocity  for  now  we  can  write 


47r 


7'''^'*<^''r^-  ^-22' 
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We  see  that  without  the  rotational  terra  the  velocity  is  a  potential  flow.  This 
fact  is  used  to  great  advantage  in  many  apphcations,  including  a  powerful  method 
of  reconstructing  the  three  dimensional  velocity  field  from  redshift  measurements, 
known  as  P0TENT.9'i° 

Finally  we  note  that  the  velocity  is  sometimes  written  in  terms  of  the  peculiar 
gravitational  acceleration,  defined  by  g  =  —V(l)/a.  We  can  then  write 

The  disregard  of  the  rotational  term  means  that  the  peculiar  velocity  and  accelera- 
tion are  parallel. 

2.2  Statistics 
Now  that  we  have  some  idea  of  how  density  fluctuations  and  peculiar  velocities 
should  act  (at  least  on  large  scales,  for  now)  we  need  to  be  able  to  compare  these 
predictions  to  observations  or  to  N-body  simulation  results.  Of  course  it  is  not  the 
density  or  velocity  fluctuations  at  any  one  point  that  are  of  primary  interest,  but 
the  statistical  distributions  of  said  quatities.  The  distributions  of  such  quantities  are 
defined  by  their  moments.  For  a  one  dimensional  variable,  say  x,  with  a  distribution 
f{x)  the  moments  are  defined  by 

(2;")-   f  dxf{x)x''.  (2.24) 

This  is  easily  generalized  to  fields,  and  the  first  few  moments  of  the  density,  p{x) , 
are  written 

(p(x))  =  po,  (2.25) 

(p(xi)p(x2))  =  Po'(l  +  ei2),  (2.26) 

{p{xi)p{x2)p{x3))  =  po^(l  +  62  +  63  +  63  +  C123).  (2.27) 
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Here  ^  is  the  irreducible  two-point  correlation  function  and  C  the  irreducible  three 
point  correlation  function.  Of  course  this  hierarchy  can  be  continued,  but  the  two 
and  three  point  functions  and  their  Fourier  transforms  are  the  most  commonly  used. 
The  two  point  correlation  function,  ^,  is  also  seen  to  be  the  second  moment  of  the 
density  perturbation  6, 

C  (P1P2)  -  PO^  ((P1-P0)(P2-P0))  , .    .  V  ,^  „„. 

02  = 5 = 5 =  \p\O2),  (2.28) 

Po^  Po^ 

where  the  subscripts  1  and  2  indicate  a  function  of  Xx  or  x^  respectively.    In  the 

second  line  we  have  used  the  fact  that  (5i)  =  (52)  =  {{p\,2)  —  Pa) / Pa  =  0,  by  equation 

[2.25]. 

The  physical  meaning  of  ^  is  easily  appreciated  when  we  look  at  its  definition 

in  terms  of  a  discrete  quantity,  such  as  the  number  counts  of  galaxies.  Assuming  no 

bias  between  mass  distribution  and  galaxy  distribution,  we  can  write  the  probability 

of  finding  a  galaxy,  which  we  assume  to  be  a  point-like  object  as  far  as  the  scales  we 

are  interested  in  are  concerned,  in  an  arbitrary  infinitesimal  volume  as  dP  =  p{x)dV. 

Then  the  probabihty  of  finding  one  galaxy  in  each  of  two  different  volumes  is 

dP  =  {p{x,)p{x2))dV,dV2  =  po'[l  +  ^{xn)]dV,dV2.  (2.29) 

We  can  see  then  that  ^,  as  its  name  suggests,  is  simply  the  correlation  between  the 
two  galaxies.  That  is,  if  ^  is  positive  for  some  X12  then  galaxies  are  clustered  at 
that  length  scale.  Given  one  galaxy  there  is  a  greater  chance  of  finding  another  at 
that  distance  than  would  be  found  by  chance,  as  determined  by  the  average  density 
of  the  universe.  Having  said  that  we  assumed  no  bias  between  mass  distribution 
and  galaxy  distribution,  we  must  at  least  acknowledge  that  such  a  bias  may  exist. 
At  large  scales  there  appears  to  be  no  bias,  so  the  two  correlation  functions  are 
equivalent,  but  on  small  scales  a  significant  bias  may  exist.  This  is  the  hkely  source 
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of  differences  between  correlation  functions  calculated  from  observations  and  those 
from  N-body  simulations  and  analytical  calculations.^^ 

It  is  often  useful  to  examine  density  fluctuations  and  their  associated  statistics 
in  Fourier  space.  In  fact,  the  second  moment  of  the  transform  of  the  density  pertur- 
bation, the  power  spectrum,  is  much  more  often  calculated  than  ^  itself.  We  define 
the  Fourier  transformed  density  fluctuation 

5ik)=   fd^xS{x)e'^'',  (2.30) 

then  the  power  spectrum  is  defined  by  ■ 

{S{ki)~S{k2))  =  [{2wfSD{ki  +  k,)]P{k),  (2.31) 

where  P{k)  is  the  power  spectrum  and  So  is  the  Dirac  delta  function.  Here  k  is 
the  magnitude  of  either  wave  vector  (as  the  Delta  function  ensures  they  are  equal.) 
Isotropy  demands  that  the  power  spectrum  is  a  function  of  magnitude  only.  With 
these  definitions,  the  power  spectrum  can  be  written  in  terms  of  the  two  point 
correlation  function  hke 

P{k)^   fd^x^{x)jo{kx),       '  (2.32) 

where  jo{kx)  =  sm{kx)/kx  is  the  zeroth  spherical  Bessel  function.  The  same  proce- 
dure can  be  undertaken  for  the  higher  moments.^' ^^ 

Finally,  we  wish  to  examine  the  moments  of  the  pecuhar  velocity  field.  As  we 
already  have  a  relationship  between  the  peculiar  velocity  and  the  density  fluctua- 
tions, it  is  straightforward  to  find  the  moments  of  the  (irrotational)  velocity  field  in 
terms  of  the  correlation  functions.  For  example, 

=     '-^  J  d'yay)/y.  (2.33) 
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gives  the  second  moment  of  the  peculiar  velocity  in  terms  of  the  two  point  correlation 
function.  Many  authors  wiU  instead  calculate  moments  of  a  variable  proportional  to 
the  velocity  divergence,  9  =  (l/d)V  •  v.  This  procedure  obviously  loses  its  utility  if 
the  rotational,  or  divergence-free,  mode  of  the  velocity  is  kept.  Many  examples  can 
be  found  in  the  hterature.^'^^'^^  Another  interesting  set  of  moments  to  calculate  are 
those  of  the  difference  between  the  velocity  of  two  galaxies.  For  example  the  first 
moment  of  one  component  of  the  density  weighted  pair  velocity  is 

where  the  subscripts  refer  to  fimctions  oi  Xi  or  aj2. 

One  reason  these  are  interesting  is  that  they  can  be  easier  to  calculate  from 
actual  data,  due  to  difficulties  in  determining  distances  to  galaxies  independent  of 
redshift.^^ 

2.3     Small-Scale  Considerations 

The  previous  perturbative  analysis  works  well  as  long  as  the  density  perturba- 
tions are  small.  {5  <  1.)  When  5  becomes  large  we  enter  the  non- linear  regime  and 
different  analytic  methods  are  required.  This  non-linear  regime  represents  objects 
which  have  collapsed,  which  we  can  recognize  as  galaxies  or  clusters  of  galaxies  with 
their  associated  dark  matter  halos.  We  typically  do  not  speak  of  individual  stars 
as  collapsed  objects  in  this  sense,  as  much  more  comphcated  physics  contributes 
at  those  scales.  There  are  two  main  ingredients  in  the  description  of  non-linear 
density  perturbations.  First  one  must  know  something  about  the  distribution  of 
sizes  of  the  collapsed  objects.  This  has  traditionally  been  accomplished  using  a 
method  known  as  the  Press- Schechter  approximation, ^''  which  has  been  improved 
upon  since  its  original  derivation.  ^^-^^  To  understand  the  Press-Schechter  approach, 
we  first  remember  that  the  density  contrast,  5,  is  a  Gaussian  random  field,  or  at  least 
is  assumed  to  be  in  this  derivation,  and  we  define  its  standard  deviation,  a  =  ((5^)^/^. 


■:r-i'^: 
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We  can  then  write  the  probabihty  that  the  density  contrast  in  some  randomly  placed 
window  exceeds  a  critical  density  for  collapse,  Sc,  as 

1        r°°  ' 

where  cr  is  a  function  of  the  window  size,  or  of  the  contained  mass,  M.  This  can  also  ; 
be  viewed  as  the  fraction  of  perturbations  within  the  window  whose  mean  density 
exceeds  Sc-  Press  and  Schechter  themselves  realized  that  this  would  only  take  into 
account  the  mass  contained  in  the  overdense  regions  of  the  universe,  and  since  the 
mass  in  the  under  dense  regions  will  eventually  accrete  onto  the  collapsed  objects, 
they  simply  multiphed  their  result  by  a  factor  of  two.  Now,  to  find  the  number 
density  of  collapsed  objects  in  the  range  dM  (actually,  it  is  typical  to  write  the 
result  in  terms  of  logarithmic  interval  d  In  M)  we  differentiate  the  above  probabihty 
and  multiply  by  the  inverse  of  the  volume,  po/M,  (which  is  simply  the  number 
density  of  collapsed  objects  of  that  size  if  all  mass  were  contained  in  objects  of  that 
size.)  This  gives  us  the  result 

dn  /2\       pod\na~^    _-i?i2 


ue-"  '\  (2.36) 


d\nM       \n J       M  dlnM 

where  u  =  bcjo.  The  value  of  the  critical  density  contrast  for  collapse  is  (5c  =  1.686. 
(See  Appendix  B  for  a  derivation  and  discussion  of  this  result.) 

Though  the  Press-Schechter  mass  function  has  been  the  standard  model  used 
since  its  inception,  it  seems  to  predict  too  many  low  mass  halos,  whereas  N-body  sim- 
ulations and  astronomical  data  suggest  that  the  low  mass  objects  are  more  strongly 
clustered  and  thus  should  properly  be  classified  as  higher  mass  objects.  There  have 
been  many  recent  efforts  to  improve  upon  the  Press-Schechter  formahsm,  including 
an  analj^ical  attempt,  based  on  the  Zeldovich  approximation^^  (using  an  ellipsoidal 
collapse  model  instead  of  the  spherical  collapse  used  by  Press-Schechter),  and  fits 
to  N-body  data.^^    Sheth  and  Tormen,^^  find  a  mass  function,  valid  for  N-body 
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simulations  of  a  variety  of  cold  dark  matter  cosmologies,  given  by 

*      -'^'^'--'s^l  +  JLUl^V^-V.,  (2.37) 


dlnM       M  d\nM        \        i/'^y  ^2^/ 

where  i/  =  a-^/^i/,^  =  0.707,  g  =  0.3,  and  A  =  0.322.  Note  that  the  values 
a  =  1,  5  =  0,  yl  =  1/2  reproduce  the  Press-Schechter  function. 

The  second  ingredient  in  describing  non-linear  perturbations  is  a  knowledge  of 
the  shape  and  nature  of  the  objects  that  have  collapsed.  This  has  been  possible 
only  through  the  analysis  of  N-body  simulations,  where  (possibly)  universal  profiles 
of  collapsed  objects  have  been  extracted.  The  first  successful  universal  profile  for 
dark  matter  halos  was  presented  by  Navarro,  Frenk,  and  White  (hereafter  NFW)  in 
1997.^^  Their  study  of  N-body  results  led  to  the  suggestion  of  a  universal  density 
profile  of  the  form 

"pT  ^  (r/r,)(l  +  r/r,)2'  ^^"^^^ 

where  pc  —  SH'^/SttG  is  the  critical  density  for  closure  of  the  universe,  Sa  is  a 
characteristic  density,  and  rg  is  a  scale  radius  which  is  different  for  each  individual 
cluster.  (The  scale  radius  is  found  by  N-body  simulations  to  be  a  rough  function  of 
cluster  mass,  a  fact  we  will  use  in  chapters  5  and  6.)  More  recent  results  of  numerical 
simulations,  performed  by  Moore  and  collaborators,  with  higher  resolution  have 
indicated  that  the  universal  profile  may  instead  be^^ 

P^  = ^^ (239) 

As  we  have  often  made  mention  of  the  results  of  N-body  simulations,  it  would 
be  appropriate  at  this  point  to  make  a  few  remarks  regarding  them.  Because  exact 
solutions  to  non-hnear  equations  of  motion  are  difficult  to  come  by,  we  are  led  to  try 
to  use  computer  simulation  to  gain  insight.  These  N-body  simulations  approximate 
the  continuous  density  field  by  a  set  of  discrete  particles.   The  simulations  simply 
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calculate  the  gravitational  force  on  each  particle  due  to  all  the  other  particles,  ac- 
celerate and  move  the  particles  appropriately  by  a  small  amount,  then  iterate  the 
procedure.  Of  course  this  procedure  quickly  becomes  intractable  when  the  number 
of  particles  becomes  large,  so  many  methods  have  been  developed  to  shorten  the 
calculations.  Typically,  as  in  the  most  popularly  used  codes,  called  particle-particle- 
particle-mesh,  or  P^M  codes,  exact  forces  are  only  calculated  for  particles  close  by, 
with  the  remaining  particles  binned  onto  a  grid.^^ 

The  use  of  discrete  particles  to  represent  an  essentially  continuous  density  cre- 
ates its  own  problems.  The  most  serious  of  these  arises  when  particles  pass  close  to 
one  another  and  scatter  at  wide  angles.  Since  each  particle  actually  represents  a  huge 
number  of  smaller  particles  spaced  far  apart,  such  collisions  should  not  take  place. 
In  order  to  correct  for  this,  forces  are  "softened"  over  appropriately  small  length 
scales.  This  effect,  along  with  the  complicated  physics  that  arises  out  of  the  colli- 
sions that  do  take  place,  such  as  shock  heating  and  radiative  cooling,  significantly 
limits  the  resolution  of  N-body  simulations. 

One  of  the  key  ingredients  in  determining  the  outcome  of  the  simulations  is  the 
selection  of  the  initial  conditions.  These  conditions  depend  on  the  cosmology  which 
one  is  trying  to  simulate,  but  typically  the  particle  positions  are  distributed  in  a 
Guassian  fashion  and  the  initial  velocities  are  projected  onto  the  growing,  radial, 
mode. 


CHAPTER  3 
VELOCITY  MODES  IN  PERTURBATION  THEORY 

3.1     First  Order  Solutions 

In  this  chapter  we  wish  to  systematically  include  rotational  velocity  terms  in 
the  structure  of  perturbation  theory.  We  begin  by  rewriting  the  equations  for  the 
evolution  of  5  and  v  from  the  last  chapter.  Equation  [2.12]  is  written  with  the  first 
order  terms  on  the  left  and  the  second  order  source  term  on  the  right, 

S+-V-v^--V-{Sv).  (3.1) 

a  a 

We  also  remember  equation  [2.17],  which  we  divide  by  the  scale  factor  a, 

where  we  have  used  the  fact  that  a  a  t^^^  and  eqns.  [2.9]  and  [2.15].  As  we  saw,  the 
first  order  solution  for  5  is 

,      6i  =  ^(x)f2/^  +  B{x)t-\  (3.3) 

where  A  and  B  are  arbitrary  functions  of  position  which  represent  the  density  per- 
turbations at  some  initial  time.  Here  and  in  the  rest  of  the  chapter  we  adopt  the 
notation  used  by  Hwang.'^"''^^  The  first  order  solution  for  v  is  obtained  by  setting 
the  right  hand  side  of  equation  [3.1]  equal  to  zero  and  putting  in  the  above  solution 
for  S.  We  thus  have 


a  a" 


=    VV" 


AAt-'/^  +  Bt-^ 
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+  Ct'^^/^  (3.4) 
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where  V  •  F  =  V  •  C  =  0.  The  last  term  is  the  rotational  mode.  For  algebraic 
ease  we  will  set  S  =  0  as  this  is  the  decaying  mode  in  the  first  order  solution  for 
5.  The  function  A  then  is  the  density  perturbation  at  some  initial  time,  A{x)  = 
S{x,to)  =  5o{x).  The  vector  C  is  the  initial  rotational  velocity,  C{x)  =  vl{x).  As 
A  and  C  describe  initial  conditions,  they  are  typically  taken  to  be  random  variables 
with  Gaussian  distribution. 

3.2     Second  Order  Solutions  in  Real  Space 
In  order  to  find  second  order  solutions  for  S  we  plug  the  first  order  solution  into 
the  right  hand  side  of  equation  [3.1], 


4  +  -V  •  ■U2  =  — V  •  (viSi). 
a  a 


(3.5) 


We  then  have 


^2  +  -V  •  V2 

a 


o 


(3.6) 


where 


El  =  V  •  (^VV-M)     and    E2  =  ~V  ■  {CA)  = -C -VA. 


(3.7) 


Doing  the  same  with  equation  [3.2]  we  have 


_5    V2        va^      _2 
dt^  a'       a3t      St"^ 


(-)  +  -^  +  ^vv-% 


V 


-VV-^^At-'^^  +  Ct-'/' 


) 


(3.8) 
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where 


Fi    =  -(VV^M)  •  V(VV-M) 

=  -iv(VV-M-VV-M)  (3.9) 

F2    =  (VV-M-V)C  +  C-V(VV-M)  (3.10) 

F3    =  -{C-V)C.  (3.11) 

Taking  the  divergence  of  equation  [3.8]  and  inserting  equation  [3.6],  we  obtain  one 
second  order  differential  equation  for  (^2,  • 

=  0,  (3.12) 

which  admits  a  solution  of  the  form 

52  =  Hr^'^^iei^  +  Jt'"\  (3.13) 

where  H,  I  and  J  are  arbitrary  functions  of  position.  Plugging  this  in  and  matching 
powers  of  t,  we  obtain  the  result 

5^  =  1{5E,  -  2V  •  Fi)t'/'  +  ^(-^2  +  V  •  F2y^^  +  (^V  •  F^)t-^/\         (3.14) 

The  first  term  is  still  the  dominant  one,  and  has  no  contribution  from  rotational 
modes.  We  will  revisit  this  term  in  a  moment,  when  we  look  at  its  Fourier  coun- 
terpart. The  second  term  describes  the  mixing  of  the  rotational  and  longitudinal 
modes.  An  important  feature  is  the  fact  that  this  mode  does  not  decay  with  time 
buts  actually  grows,  meaning  that  even  an  initially  small  rotational  component  will 
have  an  impact  on  the  future  longitudinal  mode.The  third  term  is  due  to  the  purely 
rotational  mode,  and  it  decays. 
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Putting  this  solution  back  into  equation  [3.6],  we  see  that  the  solution  for  the 
second  order  velocity  looks  like 


a 


1,      .„  8„      "   -.1/3    ,    /3  i  p   w-2/3    ,    /3  i^  X       5/3 


-(-2^1  +  -V  •  F^y^'  +  {-E^  -  -V  •  F2)r2/^  +  (-V  •  F,)t 
+     divergence-free  part.  (3.15) 

To  find  the  divergence  free  part,  we  look  at  the  curl  of  equation  [3.8]  and  note  that 

V  X  Fi  =  0.  This  implies 

where  the  superscript  t  indicates  the  transverse  (or  divergence-free)  part,  which  is 
defined  by  F^^  =  {S^  -  VaV~2V^)F„^.2i  Here  and  for  all  that  follows,  a,  /?  represent 
spatial  indices  and  thus  run  from  1  to  3.  Equation  [3.16]  admits  a  solution  of  the 
form 

^-Lr2/3  +  Mr^/^  (3.17) 

a 

where  L  and  M  are  arbitrary  functions  of  position  with  the  caveat  that  V  •  L  = 

V  •  M  =  0.  When  L  and  JVI  are  solved  for  and  inserted  into  the  total  solution  we 
have 


"^  =  vv- 


^(-2Fi  +  ^V  •  F,)t'/'  +  {\e,  -  \v  ■  F,)t-'l'  +  (^V  •  F,)r"' 
+F*t-2/^  -  3F^r^/^  (3.18) 

We  note  that  to  second  order  the  velocity  still  has  only  one  growing  mode,  which  is 
still  purely  longitudinal.  It  does,  however,  have  several  terms  constant  in  time,  one 
of  which  is  a  potential  flow  which  is  coupled  to  the  initial  rotational  mode  and  one 
which  is  a  pure  rotation,  coupled  to  the  initial  longitudinal  mode. 
3.3  Second  Order  Solutions  in  Fourier  Space 
It  is  often  convenient  to  work  in  Fourier  space.  (One  reason  for  this,  among 
many  others,  is  that  spatial  derivatives  become  products  with  the  wave- vector.)  We 
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define  the  Fourier  transform  in  the  normal  way,  i.e., 


(3.19) 
(3.20) 


To  first  order  then,  disregarding  the  decaying  density  term,  we  can  write  Si  =  A0^, 
and  j«^ 

^  =  4T^At-'/'  +  Ct-'^\  (3.21) 

To  transform  the  second  order  density,  we  look  at  each  of  the  terms  in  equation 

[3.14]  in  order.  The  first  term  contains  Ei  and  the  divergence  of  Fi,  which  can  be 

'  ?  '■  ' 

written  %  ,        • 


El    =    V  •  (^VV-M)  =  (VA)  •  (VV-M)  +  yl' 
V-Fi    =    -^V  •  V[(VV-2yl)  •  VV-M] 


(3.22) 


=    -V«V;3(V-'^)V"V^V-M  -  V^^V-MV^yl.  '   ^'  (3.23) 
The  first  term  of  equation  [3.14]  can  therefore  be  written  > 


5.„  =  t'l^ 


-A^  +  (V^)  ■  (VV-M)  +  -(V«V;3(V-'^)V"V^V-M) 


(3.24) 


which  is  a  well  known  result. ^^'^^'^^  In  Fourier  space,  products  become  convolutions, 
as  defined  by 

Fi{x) .  ..F,{x)  =  y  (^  •  •  •  ^s^oi^f'^  -  fc)A(fci)  •  --Fnik^),  (3.25) 

where  a  quantity  with  a  tilde  denotes  its  Fourier  transform  and  Sd  indicates  the 
Dirac  delta  function. 

So  we  have  ' 


~S2a  -  t^"  j ^,^,[{2.f5n{ki+k,-k)]A{ki)A{k, 


5      fci  •  fc2      2  (fci  •  ^2) 


+ 


+ 


i  fto  /  rC-i  /Co 

(3.26) 


^>,  -~ 


24 
The  second  term  contains  £'2  and  the  divergence  of  F2,  so  we  write 

E2    =    -C^VM  (3.27) 

V.F2   =   V  •  [(vv-M  •  V)C  +  C  •  V(  vv-M)] 

=    (V^C«)VaV;5V-M  +  (V,C"')(V^V"V-M)  +  C^VfsA.    (3.28) 


Thus 


and 


^2.  =  t'/' 


-CaV"^  +  (Vc,V^V-M)(V^C°  +  v°c^) 


(3.29) 


5,.  =  t^/^ 


y   (27r 


)3(2^ 


[(27r)^<5c(fei  +  fc2-fc)] 


z^C'i-M2  +  i2ii^^fci-C2 


(3.30) 


where  the  subscipts  on  A  and  C  refer  to  functions  of  fci  or  fc2.    The  third  term 
contains  only  the  divergence  of  F3,  which  is 


V  •  F.  =  -V„C^V.C° 


(3.31) 


Thus 


<52c  =  t-'/'  (  -^V.C'^V^C"  j  , 


(3.32) 


and 


4  =  r^/'  /  ^|!^t(27r)^5,,(fci  +  fc2  -  fc)]^(fci  •  C2)(fe2  •  Ci).         (3.33) 

The  final  result  for  the  Fourier  transform  of  the  second  order  density  is  of  course 
just  the  sum  of  these  terms, 


^2  =  ^2a  +  ^26  +  ^2c- 


(3.34) 


The  Foiuier  transform  of  the  second  order  velocity  is  found  in  similar  fashion. 
Looking  at  each  of  the  four  terms  in  equation  [3.18]  in  order,  we  have,  for  the  first 


^^^.:-- 
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term 


(-) 


VV 


-2 


8. 


^i/J    (_2^i  +  -V-Fi) 


7 


^    -VV-H'/^ 


8 


^^2  +  ^(Vyl)-(VV-M) 


+  -V«V^(V-M)V°V^(V-2A) 


(3.35) 


Its  Fourier  transform  is  then 

'V2 


^).  =  -I'-VS^^K^"''^""""*-^"^'^^ 


X 


2      2  fci  •  fca        8  (fei-fca)^ 
7^3     iti      + 


(3.36) 


1,2  ^i-        "^1^2 

This  is  the  full  solution  when  the  rotational  mode  is  ignored  and  is  well  known  in 
the  hterature.^'^'^'^^  The  second  term  is 


\  a  /b 


t-'^HlE,-\w.F,) 


whose  Fourier  transform  is 


V2 


,  (3.37) 


=   -'"^'F/^il|K^^''*''('='^^-^)l 


X 


z2Ci  •  k2A2  +  iAj^^^ki  ■  C2 


(3.38) 


The  third  term  is 


(3.39) 
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whose  transform  is 


V2 


;(«!,■  C2)(fc2-Ci) 


The  final  term  is  the  divergence  free  part, 


p)  =  F*r2/3  -  3F*r^/^ 


(3.40) 


(3.41) 


Using  the  previous  definitions,  we  can  write  this  term  in  component  notation  as 
\  a  /  da  "■ 

-3r^/^[C^V^Ca  -  VaV-^V^iC-yV^Cp)].  (3.42) 

The  Fourier  transform  of  this  term  is 


V2 

a 


=  *-"iSts$^si(^^no(M^k.-k)] 


X 


«=' • '=^6.  +  i.^i^fc,  -  * .  ( ^c,  +  i,^i^fc. 


,2       —  z    .    — ^        ,2        '"^        '■'      I         i.2 
1  2  \  1 

(Ci  •  fc2)C'2  -  ^  (Ci  •  ^2)02 


1^2 


(3.43) 


Written  in  this  form  it  is  easy  to  see  that  this  contribution  is  purely  perpendicular 
to  the  wave-vector  k. 

3.4     Third  Order  Solutions  in  Real  Space 
Solutions  to  third  order  can  be  obtained  in  much  the  same  way  as  second  order 
solutions.  They  are  necessary  in  order  to  calculate  quantities  which  are  sometimes 
of  interest,  such  as  the  four-point  function  or  its  Fourier  analog,  the  trispectrum.^^'^^ 
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To  third  order,  equation  [3.1]  becomes 
a  \  a     J  \  a     / 


=    -V 


'\v-F^t--"^  +  ^-{-E,  +  V-F^)t"^ 


-V  •  [(ylt2/3)|vv-2(fV3l(_2£;j  +  ^y  •  Fj) 

+r2/3(3^^  -  i V  •  F2)  +  t-"^{\v  ■  F3))  +  F\r^'^  -  3F*r^/^}l 


^Po  +  Pi  +  r^Pa  +  r'Pa, 


(3.44) 


where 


Po    = 


-^■\ 


{^E^-2V-Fi)VV -''{-- A) 


8. 


+^VV-^(-2£;i  +  -V  •  Fi) 

.    .  .  O 


(3.45) 


Pi    = 


{-E^  +  V  •  F2)  vv-2(-|a)  +  Avv-2(^£;2  -  ^v  •  F2) 


1 


+-(5^1  -  2V  •  Fi)C 


-V 


(^V  •  F3)VV-2(-^^)  +  ^(-^2  +  V  •  F2)C 


+^VV-2(-V  •  F3)  -  3^F^ 


-V 


(^V  •  F3)C 


(3.46) 


(3.47) 
(3.48) 


Po  contains  all  pure  longitudinal  couphngs,  Pi  and  P2  contain  mixed  couplings,  and 
P3  contains  the  purely  rotational  coupling. 

Similarly,  to  third  order,  equation  [3.2]  becomes 


dt    a 


a  3 


a 


V— 
a        a         a 


=    Qo  +  Qit~'  +  Q2t~^  +  Q3t~^   (3.49) 
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where 


Qo    =    -[VV-^(-^yl).V][VV-^^(-2£;i  +  ^V-F0]  '     r 

-[(vv-^(^£;2  -  ^v  •  F,)  +  F\) .  v][vv-\-Ia)] 

-(C-V)[VV-2^(-2£;i  +  ^V.Fi)] 
-[(VV-^^[-2Ei  +  ^V.Fi])-V]C, 

-[(VV-^[^V  •  F3]  -  3F^) .  V][VV-^(-^^)] 

-[(vv-2[^£;2-^v.F2]  +  F^).v]c, 

Q3    =    -(C-V)(-3F|)-(-3F*-V)a  (3.50) 

As  before,  we  combine  the  divergence  of  equation  [3.44]  with  equation  [3.49]  to 
produce  one  second  order  differential  equation  for  ^3.  This  is 

k  +  \t~H^  -  y-%  =  (-^Po  +  V  •  Qo)  +  (^Pi  +  V  •  Q,)r' 
H-\P2  +  V  •  Q^)t-^  +  (^P3  +  V  •  Q^)t-\  (3.51) 

which  admits  a  solution  of  the  form 

5^^Re^St  +  U  +  Vt-\  (3.52) 

Plugging  in  this  solution  and  matching  powers  of  t  gives  us  the  final  solution 

-\i-\P2  +  V  •  Q2)  -  ^(^^3  +  V  •  Q,)t-\  (3.53) 
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The  first  term  is  the  purely  longitudinal  mode,  and  is  again  the  fastest  growing  term. 
The  second  term,  however,  is  also  growing,  and  contains  a  mixture  of  longitudinal 
and  rotational  couplings.  The  third  term  is  a  constant  in  time  and  is  also  a  mixture  of 
couplings.  The  fourth,  and  still  decaying,  term  contains  purely  rotational  couplings. 
Putting  this  solution  into  equation  [3.44],  we  find 


a 


ti^Po  +  V  •  Qo)  -  (Pi  +  V  •  Qi)  +  r'P2  +  t-^lPs  +  V  •  Q3) 
b  o 


+     divergence- free  part. 


(3.54) 


To  find  the  divergence-free  part,  we  note  that  V  x  Qq  =  0  and  look  at  the  purely 
rotational  part  of  equation  [3.49], 


at  \  a  I       a  6 


which  admits  a  solution  of  the  form 


V 


-A  =  W  +  Yr^  -h  zt- 


(3.55) 


(3.56) 


Matching  powers  of  f ,  we  obtain  the  final  result 
H"7"^o-t-  V 

D 

3 


=    VV 


tC^Po  +  V  •  Qo)  -  (Pi  +  V  •  Qi)  +  r'p^  +  t-\\P:,  +  V  •  Q3) 
6  3 


+^Qi  +  3Q^r  1 


-^Qlt-\ 


(3.57) 


The  fastest  growing  mode  of  this  solution  is  the  purely  longitudinal  peice,  but  there 
are  several  other  growing  terms  as  well  (remembering  that  the  a  in  the  denominator 
of  the  left  side  is  proportional  to  t^/^.)  The  second  term  is  a  potential  flow,  grow- 
ing in  time,  which  derives  from  a  mixture  of  initially  longitudinal  and  rotational 
modes.  The  other  growing  term,  (3/4)Q*,  is  a  purely  rotational  mode  deriving  from 
a  mixture  of  initially  longitudinal  and  rotational  modes.  Thus  it  is  at  third  order  in 
perturbation  theory  that  we  find  growing  terms  in  the  peculiar  velocity  associated 
with  the  initial  rotational  flow.  This  should  be  somewhat  important  in  quasi-linear 
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regimes  where  5  is  approaching  unity,  although  perturbation  theory  at  third  order 
and  higher  has  a  fairly  small  window  of  applicability. 

3.5     Evolution  of  Vorticity 
As  a  corollary  to  the  discussion  above,  we  look  at  the  evolution  of  vorticity 
in  perturbation  theory.  We  define  the  vorticity,  C,  to  be  the  curl  of  the  transverse 
velocity.  Remembering  the  first  order  solution  for  the  pecuhar  velocity, 


V 


-V  fd'x'Six')-^  +  ^  =  «^^o  +  -,  (3.58) 

J  \x'  —  x\  a  a 


in 
we  find  that  the  first  oder  vorticity  is 


a        a 


(3.59) 


Thus  to  first  order,  the  vorticity  decays  as  the  inverse  of  the  scale  factor,  a.  To 
explore  the  time  dependence  of  the  next  order  term,  we  start  with  equation  [2.17] 
from  the  last  chapter.  We  then  use  the  vector  identity 

V  X  (V  X  v)  =  V(-u2)  -  {v  ■  V)t;,  (3.60) 

which  gives  us 

Id  ,     ^      1„  /I   ,\       1         ,  1 

Taking  the  curl  of  this  equation  and  writing  it  in  terms  of  C,  we  have 

^(aC)  =  V  X  (v  X  C).  (3.62) 

Using  the  first  order  velocity  and  vorticity  as  the  source  for  the  second  order  vorticity, 
we  find 

^{a^^)  =  V  X  (t,(i)  X  e^)  =  DVx  (770  x  Co)  +  ^V  x  (vot  x  Co)-        (3.63) 
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The  first  term  dominates  for  late  times  and  so,  disregarding  the  second  term,  we  see 
that 

C^^'  =  ^V  X  (770  X  Co)-  (3.64) 

The  interesting  aspect  of  this  result  is  that  it  is  the  interaction  between  the  first 
order  vorticity  and  the  first  order  longitudinal  mode.  This  of  course  just  echoes  the 
result  of  equation  [3.14]  and  the  fact  that  as  an  object  collapses  any  small  rotational 
velocity  is  amphfied  due  to  angular  momentum  conservation. 


CHAPTER  4 
THE  CMB  AND  THE  SUNYAEV-ZELDOVICH  EFFECT 

4.1     Introduction 

One  possible  place  in  which  the  effect  of  rotational  velocities  may  be  important 
is  in  the  scattering  of  the  cosmic  microwave  background  (CMB)  by  foreground  ion- 
ized gas.  One  of  the  major  predictions  of  the  now  almost  universally  accepted  hot  big 
bang  theory  is  the  existence  of  background  relic  radiation.  This  is  the  radiation  seen 
from  the  surface  of  last  scattering.  Early  in  the  universe  (z  >~  1000)  temperatures 
were  high  enough  that  neutral  atoms  could  not  form.  At  this  time  photons  scattered 
readily  off  the  densely  packed  charged  particles.  At  about  z  =  1000  the  tempera- 
ture dropped  enough  to  allow  the  production  of  neutral  hydrogen  and  photons  were 
"frozen  out."  It  is  from  this  epoch  that  we  receive  the  background  radiation.  The 
existence  of  this  radiation  was  first  predicted  by  Gamow,  Alpher,  and  Herman  in  he 
1950's,  and  the  search  for  it  began  in  earnest  in  the  mid  1960's.  It  was  discovered, 
rather  accidentally,  by  Arno  Penzias  and  Robert  Wilson  in  1965.^^  Since  that  time 
a  great  many  measurements  of  the  temperature  and  polarization  of  the  CMB  have 
been  made.  The  intensity  of  the  CMB  has  been  measured  very  accurately  over  more 
than  three  orders  of  magnitude  in  frequency  and  the  spectrum  is  found  to  be  that 
of  a  blackbody  with  a  temperature  of  2.725  ±  0.002K.^*^  Figure  [4-1]  shows  how 
well  the  data  fit  a  Planckian  spectrum.  One  very  intriguing  aspect  of  the  CMB  is 
its  near  total  isotropy.  The  temperature  in  every  direction  is  identical  to  a  part  in 
10~^.  This  thermal  equihbrium  over  so  many  causally  disconnected  regions  is  now 
reasonably  well  explained  by  theories  of  inflation.  (In  fact,  it  was  one  of  the  principle 
motivations  for  the  development  of  inflationary  theories.) 
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Wavelength  (cm) 
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1000 


Figure  4-1.  Measurements  of  the  intensity  of  the  CMBR  as  a  function  of 
frequency.  The  soHd  hne  represents  a  perfect  blackbody  at  2.73 
K.  The  figure  is  taken  from  Smoot.^^  (See  references  therein  for 
sources  of  data.) 
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4.2     Anisotropies  in  the  CMB 

Perhaps  more  interesting  than  the  problem  of  explaining  the  isotropy  of  the 
CMB  is  explaining  its  anisotropies.  These  anisotropies  are  tjrpically  expressed  in 
terms  of  an  expansion  in  spherical  harmonics, 

AT" 

—  {e,^)  =  Y,CilrnYlm{OA).  (4.1) 

Im 

We  can  then  define  an  angular  correlation  function  by 

C{0)  =  (^(n0^(n2))  =  Yl  ^^'^'(^'i  •  "'2),  (4.2) 

where  Pi  are  the  Legendre  poljoiomials  and  Q  =  (|aj"p),  is  the  angular  power 
spectrum.  The  first,  and  dominant,  anisotropy  is  the  /  =  1  dipole  moment.  This 
is  simply  interpreted  as  due  to  the  motion  of  the  earth  relative  to  the  "cosmic  rest 
frame"  as  defined  by  the  CMB.  Utilizing  a  simple  special  relativistic  calculation, 
Peebles  and  Wilkinson, ^^  first  noted  that  the  temperature  measured  by  an  observer 
in  a  frame  moving  with  respect  to  the  homogeneous  background,  at  an  angle  9  with 
respect  to  the  motion  is 

.,2, 


^w^StS'  (^•^) 


which  can  be  expanded  as 


-cos^+-- 
c  2c^ 


T'{9)  =  T[1  +  -  cos^  +  -^  cos  20  +  0{v%  (4.4) 


The  amplitude  of  the  dipole  anisotropy,  as  measured  by  the  Cosmic  Background 
Explorer  (COBE)  satellite,  is  AT/T  =  3.372  ±  O.OUmK.'^^  This  implies  a  solar 
system  velocity  of  371  ±  Ikm/s  towards  galactic  coordinates  {£,b)  =  (264°.  14  ± 
0.30, 48°.26±  0.30).  '.; 

4.3     Primary  Anisotropies  '    " 

After  the  dipole  anisotropy,  which,  as  we  saw,  is  at  the  AT/T  ^  10~^  level, 
is  accounted  for,  the  next  level  of  anisotropies  are  found  at  AT/T  ^  10"^.  These 
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anisotropics  are  believed  to  be  due  to  inhomogeneities  in  the  surface  of  last  scat- 
tering. The  effect  of  density  perturbations  present  at  that  time  can  be  divided  into 
three  categories.  The  first,  known  as  the  Sachs- Wolfe  effect,  is  due  to  the  photons 
climbing  out  of  the  gravitational  potentials  of  the  density  perturbations.  The  mag- 
nitude of  the  effect  is  really  the  combination  of  two  effects.  First,  as  the  photons 
chmb  out  of  the  well,  they  are  redshifted,  or  cooled,  which  gives 

AT      5(/) 

The  photons  are  also  heated,  due  to  the  associated  time  dilation,  and  you  seem  to 
be  looking  at  a  younger,  hotter  universe.  This  effect  is 

At      5(j)         AT         2(50 
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since  T  oc  1/a  and  a  oc  t^^^.  Thus  the  total  Sachs- Wolfe  effect  is  given  by 

AT       (50 


(4.6) 
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The  second  effect  related  to  density  perturbations  at  the  surface  of  last  scattering 
is  called  the  adiabatic,  or  intrinsic,  perturbation  effect.  It  is  simply  due  to  the 
fact  that  denser  regions  would  recombine  later,  thus  look  hotter.  This  effect  yields 
a  temperature  perturbation  equal  to  the  density  perturbation  responsible  for  the 
late  recombination.  The  final  important  source  of  primary  anisotropy  is  simply  the 
doppler  shift  of  the  photons  coming  off  the  moving  plasma  at  last  scattering.  This 
yields  a  temperature  perturbation  of 

ST      Sv  -f  ,     , 

where  f  represents  the  unit  vector  along  the  hue  of  sight.  The  combination  of  these 
effects  lead  to  the  famous  acoustic  peaks  in  the  angular  power  spectrum  of  the 
CMB  anisotropics.    Before  recombination,  adiabatically  ovcrdense  regions  initially 
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collapse  due  to  gravity,  but  are  then  pushed  apart  by  radiation  pressure,  setting  up 
acoustic  waves.  These  standing  waves  are  frozen  out  at  decoupling  and  are  seen,  via 
the  aforementioned  effects,  in  the  angular  power  spectrum.  The  first,  and  largest, 
peak  is  associated  with  a  length  scale  equal  to  the  Hubble  radius  at  the  time  of 
recombination.  The  size,  shape,  and  position  of  this  peak,  as  well  as  the  subsequent 
peaks,  teU  us  a  great  deal  about  allowed  cosmologies  ajid  the  values  of  cosmological 
parameters.  For  example,  as  previously  mentioned  in  chapter  1,  the  location  of  the 
first  peak  at  £  ^  200  strongly  indicates  that  the  universe  is  flat,  or  very  close  to  it. 

As  mentioned,  the  first  significant  detection  of  the  CMB  anisotropics  was  by 
the  COBE  satellite.  Since  that  time,  many  ground  and  balloon  based  measurements 
have  been  attempted.  In  2003  the  most  significant  of  these  missions,  the  Wilkinson 
Microwave  Anisotropy  Probe  (WMAP)  completed  its  first  year  of  data  acquisition. 
Figure  [4-2]  shows  the  angular  power  spectrum  measured  by  WMAP.  The  first 
acoustic  peaJj  is  very  clearly  seen.  The  second  peak  is  also  clear  and  the  beginnings  of 
a  third  are  hinted  at.  Figure  [4-3]  shows  the  results  of  a  number  of  other  experiments 
and  how  they  relate  to  the  WMAP  data.  The  agreement  among  the  variety  of 
diff^erent  measurements  shows  the  robustness  of  the  results.  WMAP  provides  a  very 
accurate  picture  of  the  anisotropy  spectrum  from  about  £  =  10  to  about  (.  =  800. 
The  large  error  bars  at  very  low  (.  are  due  to  the  fact  that  we  have  only  one  universe 
to  observe,  and  at  large  angular  scales  one  cannot  average  over  significantly  many 
independent  pieces  of  sky.  (This  is  equivalent  to  saying  that  there  are  only  2£+ 1  m's 
to  average  over.)  The  fimits  at  high  ^'s  are  due  to  resolution  hmits  for  the  WMAP 
sateUite.  The  European  Space  Agency  plans  launch  of  the  PLANCK  satellite  in  or 
around  2007,  which  hopes  to  extend  measurements  to  perhaps  £  ~  1200. 

4.4     Secondary  Anisotropies 

The  primary  temperature  anisotropies,  due  to  the  nature  of  the  surface  of  last 
scattering,  dominate  the  angular  power  spectrum  at  small  and  medium  £'s,  but 
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Figure  4-2.  The  points  represent  the  data  obtained  from  the  first  year  run 
of  the  WMAP  satellite,  along  with  \o  error  bars.  The  line  rep- 
resents a  best  fit  ACDM  model,  and  the  gray  shading  represents 
the  \a  error  due  to  cosmic  variance.  Taken  from  Hinshaw  et 
al.29 
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Figure  4-3.  A  collection  of  recent  data  from  various  CMB  anisotropy  mea- 
surement experiments.  The  line  is  again  the  best  fit  ACDM 
model  based  on  the  first  year  WMAP  results.  Taken  from  Hin- 
shaw  et.  al.^^ 
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they  are  by  no  means  the  only  contributors  to  the  spectrum.  At  smaller  scales, 
or  larger  ^'s,  secondary  anisotropies  are  expected  to  dominate  the  angular  power 
spectrum.  These  effects  include  the  Integrated  Sachs- Wolfe  (ISW)  effect,  which, 
although  it  is  a  secondary  effect,  is  expected  to  only  be  seen  at  large  angular  scales, 
and  the  thermal  and  kinetic  Sunyaev-Zeldovich  (SZ)  effects,^°'^^  which  should  be 
the  major  contributions  to  the  spectrum  at  ^'s  greater  than  about  1200.  The  ISW 
effect  is  the  secondary  analog  to  the  primary  Sachs- Wolfe  effect.  When  a  CMB 
photon  falls  into  a  potential  well  on  its  journey  to  us  it  is  heated,  or  blueshifted, 
then  cooled,  or  redshifted,  as  it  climbs  back  out.  If  the  depth  of  the  well  remains 
constant  while  the  photon  traverses  it,  then  the  redshift  exactly  cancels  the  blueshift 
and  no  trace  is  seen  in  the  CMB.  However,  if  the  depth  of  the  well  changes  over 
the  time  that  the  photon  takes  to  traverse  the  well,  then  there  is  a  net  change  in 
the  photon  temperature.  The  photon  may  traverse  many  such  regions  on  its  path 
to  us,  thus  the  word  "integrated"  in  its  name.  The  ISW  can  further  be  split  into 
two  contributions,  an  "early"  effect,  due  the  changeover  from  a  radiation  dominated 
imiverse  to  a  matter  dominated  universe,  and  a  "late"  effect,  due  to  the  changeover 
from  a  matter  dominated  universe  to  a  A  dominated  one.  (Thought  to  happening 
around  now.)  The  origin  of  the  early  ISW  lies  in  the  fact  that  when  radiation 
is  controlling  the  expansion  of  the  universe,  potentials  on  the  scale  of  the  sound 
horizon  are  damped,  giving  photons  crossing  them  an  overall  boost  in  energy.  The 
early  effect  provides  a  contribution  to  the  angular  power  spectrum  in  the  region  of 
the  first  acoustic  peak,  at  low  £'s.  This  effect  has  been  seen  by  comparing  data 
from  WMAP  to  data  from  the  Sloan  Digital  Sky  Survey  (SDSS).^^  The  late  effect  is 
expected  to  be  strongly  damped  at  medium  and  higher  ^'s  as  well,  as  the  photons 
pass  through  a  great  many  regions  of  overdensities  and  underdensities.  It  is  hoped 
that  detections  of  the  late  ISW  effect  at  low  Cs  will  provide  more  information  about 
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the  dark  energy  content  of  the  universe,  as  it  is  sensitive  to  the  time  of  switch-over 
from  matter  domination  to  A  domination.  :, 

4.5  Thermal  and  Non-Thermal  Sunyaev-Zeldovich  Effects 
At  small  angular  scales,  corresponding  to  ts  over  about  1200,  the  thermal  and 
kinetic  SZ  effects  should  be  dominant.  The  SZ  effect  is  the  result  of  the  interaction 
of  the  CMB  photons  with  regions  of  charged  particles  between  the  observer  and  the 
surface  of  last  scattering.  The  thermal  SZ  effect  is  an  inverse  Compton  scattering  of 
the  CMB  off  thermal  electron  populations.  The  result  is  a  temperature  fluctuation, 
which,  at  low  frequencies,  hke  those  of  most  of  the  CMB  photons,  is 

AT  f  ksT, 

— ^  =  -2  /  -Ueardl.  (4.9) 

Recent  analysis  of  the  data  obtained  by  the  WMAP  satellite  confirms  that  the 
SZ  effect  is  swamped  by  the  primary  anisotropics  for  all  angular  scales  available  to 
it,  corresponding  to  ts  up  to  roughly  900.^^  Though  there  is  not  yet  a  large  amount 
of  angular  power  spectrum  data  taken  at  high  £'s  (see  figure  [4-3]),  at  least  one 
mission,  the  Cosmic  Background  Imager,  has  measured  power  up  to  i  ^  3500  which 
seems  to  be  consistent  with  the  power  expected  from  the  thermal  SZ  eflFect.^-^^  More 
accurate  measures  of  this  power  may  be  useful  in  constraining  the  values  of  some 
cosmological  parameters.^^ 

There  is  an  additional  effect  due  to  scattering  off  non-thermal  electron  popu- 
lations, such  as  would  be  found  in  cluster  radio  halo  sources.  Computation  of  this 
effect  is  made  difficult  by  the  highly  relativistic  nature  of  the  electrons,  but  the 
low  density  of  such  electrons,  and  the  relative  rarity  of  such  populations,  lead  to 
insignificant  contributions  to  the  angular  power  spectrum.  (Though  the  search  for 
the  non-thermal  SZ  effect  in  specific  clusters  can  prove  useful.) 
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4.6     Kinetic  Sunyaev-Zeldovich  Effect 

There  is  an  additional  temperature  fluctuation  induced  in  the  CMB  if  the  fore- 
ground material  is  in  bulk  motion.  This  effect,  called  the  kinetic  SZ  effect,  is  simply 
due  to  the  doppler  shift  induced  by  the  motion.  The  temperature  fluctuation  along 
a  line-of-sight  direction,  7,  is  given  by 

AT        f  V 

—  =       die   ^(TTUel  ■  -,  (4.10) 

where  r  is  the  optical  depth  to  scattering  through  the  electron  cloud,  n^  is  the 
number  density  of  electrons  in  the  cloud,  and  v  is  the  pecuUar  velocity  of  the  cloud. 
Measurement  of  this  effect  for  particular  clusters  can  be  used  to  calculate  the  ra- 
dial component  of  its  peculiar  velocity,  which  is  of  significant  interest,  as  the  method 
does  not  require  an  independent  distance  measurement.^^  (It  does,  however,  require 
an  accurate  estimation  of  optical  depth  and  electron  density.)  In  addition,  one  must 
disentangle  it  from  the  thermal  effect,  as  well  as  other  small-scale  anisotropics.  This 
is  typically  done  by  making  use  of  the  fact  that  the  kinetic  SZ  eflFect  has  a  different 
frequency  dependence  than  the  thermal  effect. ^^ 


CHAPTER  5 
KINETIC  SZ  ANGULAR  POWER  SPECTRUM  INGREDIENTS 

5.1     Introduction 

In  this  chapter  we  lay  out  the  ingredients  for  calculating  the  angular  power 
spectrum  of  CMB  anisotropics  due  to  the  kinetic  Sunyaev-Zeldovich  effect.  The 
angular  power  spectrum  can  be  written  in  terms  of  a  projection  of  the  actual  three- 
dimensional  power  spectrum,  which  we  will  explore  in  the  first  few  sections  of  the 
chapter.  We  will  compare  several  of  the  popular  models  used  to  build  these  power 
spectra.  The  final  sections  will  include  calculations  of  the  angular  spectrum  in  both 
the  perturbative  regime  and  the  nonhnear  regime,  without  inclusion  of  the  rota- 
tional modes.  This  gives  us  the  necessary  background  for  the  calculation  including 
rotational  modes  in  the  next  chapter. 

5.2     Linear  Power  Spectrum 

One  of  the  main  ingredients  in  all  the  calculations  to  follow  is  a  hnear  order 
density-density  power  spectrum.  The  most  commonly  used  model  is  one  originated 
by  Bardeen  et  al.^^  The  authors  matched  a  fitting  function  to  results  from  numer- 
ical simulation  to  produce  an  analytic  form  for  the  spectrum.  A  slightly  improved 
version,  used  by  Ma,^'^  and  to  be  used  as  the  input  linear  power  spectrum  in  all  that 
follows  is  given  by 
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When  the  cosmological  parameters  ^-nx-i^K.^h  appear  without  functional  de- 
pendence on  a,  they  refer  to  the  values  today.  Which  values  to  use  are  determined 
by  the  cosmology  one  wants  to  explore.  The  values  used  in  the  calculations  that 
follow  are  based  on  a  ACDM  cosmological  model  and  values  found  by  recent  ob- 
servations,'' and  they  can  be  found  in  Appendix  C.  The  values  of  the  numerical 
constants  are  given  in  table  [5-1]. 


Table  5-1.  Linear  Power  Spectrum  Parameters 


Parameter 

Value 

OL\ 

2.34 

0C2 

3.89 

OLZ 

16.1 

a^ 

5.46 

"5 

6.71 

A 

0.785 

/?2 

0.05 

^3 

0.95 

A 

0.169 
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The  factor  of  /c"  reflects  the  assumption  that  the  initial  power  spectrum,  (at 
decouphng)  was  a  scale-free  Guassian  process  and  thus  the  initial  spectrum  was  a 
simple  power  law.  We  typically  define  a  dimensionless  version  of  the  power  spectrum, 
given  by,  A  =  k^P{k)/2iT'^.  Figure  [5-1]  shows  a  plot  of  the  dimensionless  linear 
power  spectrum  given  above. 

5.3     Nonlinesir  Power  Spectrum 

Here  we  present  a  model  for  the  nonlinear  power  spectrum  first  constructed 
by  Ma  and  Fry."*^  Before  moving  on,  we  note  that  there  are  other  persciptions  for 
describing  the  nonhnear  power  spectrum,  the  most  notable  due  to  Zhang  et  al.,^ 
but  we  follow  the  Ma  and  Fry  description. 

As  we  noted  in  chapter  2,  the  two  main  ingredients  in  this  construction  are  the 
halo  mass  function  and  the  halo  density  profile.  We  use  three  different  models  for 
the  halo  density  profile.  The  first  two  axe  the  NFW  and  Moore  profiles  mentioned 
in  chapter  2.  Following  the  convention  of  Ma  and  Fry  we  write  the  profiles,  given 
by  equations  [2.38]  and  [2.39],  as 

^  =  5au{r/r,),  (5.11) 

Po 
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A  commonly  introduced  parameter,  which  we  can  use  to  write  da  and  rg  in  terms  of 
contained  mass,  is  the  concentration  parameter,  c  =  r2m/fsi  where  r2oo  is  the  radius 
within  which  the  density  is  200  times  the  mean.  This  radius  can  be  written  in  terms 
of  mass  by,  M  —  SOOTrpo^'ioo/^-  Thus  we  can  write  the  scale  radius  as 


=  ir. 


3M 


c  VSOOttpo 


1/3 


(5.13) 
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Figure  5  1.  A  plot  of  the  dimensionless  linear  power  spectrum,  A/  = 
k^P{k)/2Tr^,  for  use  in  upcoming  calculations.  This  represents 
the  linear  power  spectrum  today,  a  =  1. 


m 

The  characteristic  density,  5a  can  also  be  written  in  terms  of  the  concentration 
parameter,  and  is  given  by 

X  200c3 

Oa,NFW      —      „ri     /,     , ^ 777— TT  (5.14) 

3[ln(l  +  c)  -  c/(l  +  c)\  ^        ^ 

'^"•^    =    ln(l  +  c3/2)-  (^-15) 

For  the  plots  in  this  section  we  use  a  model  for  the  concentration  parameter  based 
on  N-body  simulation  data  of  Jing  and  Suto"*^  given  by 


c=  < 


5a{M)     :     (M/Mq)  <  W^ 

(5.16) 
9a{M)     :     (M/Mq)  >  10^^ 


Mq  =  1.989  X  10^°  kg  is  one  solar  mass.  The  third  halo  profile  we  will  use  is 
simply  a  combination  of  the  previous  two,  as  suggested  by  Jing  and  Suto.'*^  Their 
simluations  imply  that  the  NFW  profile  is  appropriate  for  cluster  halos  and  the 
Moore  profile  is  appropriate  for  galactic  halos  and  the  halos  of  small  groups  of 
galaxies.  The  cutoff  between  the  mass  of  a  cluster  and  that  of  a  small  group  of 
galaxies  is  about  M/Mq  =  10^'*,  so  the  halo  profile  we  refer  to  as  "mixed"  is  Moore's 
for  (M/Mq)  <  10^^  and  NFW's  for  (M/Mq)  >  lO^l 

The  other  important  ingredient  is  the  halo  mass  function.  In  this  section,  we 
will  look  at  both  the  Press-Schechter  mass  function,  given  in  equation  [2.36]  and  the 
Sheth-Tormen  mass  function,  given  in  equation  [2.37].  The  principle  ingredient  in 
each  of  the  mass  functions,  as  we  saw  in  chapter  2,  is  a,  the  standard  deviation  of 
the  density  contrast,  which  is  a  function  of  the  radius  of  the  sphere  within  which  we 
are  interested.  Written  in  terms  of  the  power  spectrum,  we  have 

-^(^)  =  /     ^^^P{k)W\kR),  (5.17) 
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where  W  is  the  Fourier  transform  of  the  top-hat  window  function, 

-.,  ,       3(sin  X  —  a:  cos  x) 
W{x)  =  -^ '-.  (5.18) 

The  top  hat  window  function  has  the  value  one  within  the  radius  R  and  zero  every- 
where else  and  is  used  to  remove  the  need  to  worry  about  tiny  fluctuations  outside 
the  region  of  interest.  This  can  also  be  expressed  as  a  function  of  enclosed  mass 
using  the  relation,  M  =  4TrpoR^/3.  Figure  [5-2]  shows  a  plot  of  a  vs.  M  for  the 
input  hnear  power  spectrum  given  in  equation  [5.1]. 

Finally,  we  must  put  this  together  to  make  a  power  spectrum.  We  first  note 
that  we  can  write  the  spectrum  as  a  sum  of  two  distinct  contributions,  a  one-halo 
term,  from  particles  within  the  same  halo,  and  a  halo-halo  term,  from  particles  in 
two  separate  halos.  From  Pry  and  Ma,^^  we  can  write  these  terms  as 

Pi,(A;)    =    j  dM~[r%u{kr,)]\  (5.19) 


n2 


P2h{k)    =     \j  dM^rl5au{krs)h{M) 
where  b{M)  is  a  bias  factor  given  by  Jing  (1998),^^*  as 

2  _  1  \     /    1  \  0.06-0.02n 


Pb„(A),  (5.20) 


where,  as  before,  u  =  Sc/a{M).  This  bias  factor  is  introduced  to  account  for  the 
fact  that  dark  matter  halos  cluster  differently  than  the  general  mass  density  field. 
Clusters  form  at  the  peaks  of  the  density  fields,  thus  the  bias  accounts  for  the 
difference  between  the  statistics  of  the  the  peaks  and  those  of  the  general  field. 
Figure  [5-3]  shows  the  behavior  of  b{M).  Its  value  is  near  unity  for  masses  of  the 
most  common  clusters.  It  becomes  very  large  for  very  large  mass  clusters,  but  these 
are  exceedingly  rare  occurances.  (Which  is,  of  course,  the  very  reason  the  bias  is  so 
high.)  The  second  factor  in  the  bias  parameter  is  a  phenomenological  addition  to  the 
original  formula  based  on  N-body  results  and  the  Sheth-Tormen  perscription,  and 
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Figure  5-2.  The  linear  rms  fluctuation  of  the  density  contrast  as  a  function 
of  the  mass  within  a  top  hat  window.  The  mass  is  expressed  in 
multiples  of  one  solar  mass  (IM©  =  1.989  x  10^°%). 
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it  is  therefore  neglected  when  the  Press-Schechter  formaUsm  is  used  in  the  figures 
below. 

Figures  [5-4],  [5-5],  and  [5-6]  show  the  nonlinear  power  spectra  for  the  NFW, 
Moore,  and  mixed  halo  profiles,  respectively.  In  each  plot  the  spectrum  is  plotted 
twice,  once  for  the  Press-Schechter  mass  function  and  once  for  the  Sheth-Tormen 
mass  function.  Also,  the  linear  power  spectrum  of  figure  [5-1]  is  included  for  com- 
parison. 

Figure  [5-7]  shows  the  nonlinear  power  spectrum,  using  the  mixed  halo  profile 
and  ST  mass  function,  as  we  will  henceforth,  broken  down  into  its  one-halo  and  halo- 
halo  contributions.  Of  course  the  halo- halo  term  dominates  at  small  k,  corresponding 
to  large  distances,  or  large  masses,  and  the  one-halo  term  dominates  at  large  k, 
corresponding  to  small  scales. 

5.4     Angular  Power  Spectrum  Machinery 

Now  that  we  have  three  dimensional  power  spectra  for  both  Unear  and  nonlinear 
regimes,  we  need  to  know  how  to  use  these  to  calculate  the  angular  spectrum  of 
CMB  temperature  fluctuations  due  to  the  kinetic  SZ  effect.  This  machinery  is  due 
to  Ma  and  Fry.^^  First  we  remember  from  chapter  4  that  the  fractional  temperature 
anisotropy  due  to  the  kSZ  effect  is 

AT 


"(7) 


=   /  dle-^neOTl  ■  {v/c).  (5.22) 


T 

We  write  the  electron  density  as  ng  =  ne{l  -\-  5)  where  the  mean  electron  density 
can  be  written,  fie  =  Xe^bPci^  +  zY/nip.  Here  Xg  is  the  ionization  fraction,  rup  is 
the  mass  of  the  proton,  and  pc  =  SH^/SttG  is  the  critical  density  for  closure  of  the 
universe.  We  can  then  write  the  temperature  fluctuation  as 

AT  f 

—  {j)  =  J  die  ^  near  J  ■  (q/c),  (5.23) 


50 


106 


10^ 


10* 


S   1000 


100 


10 


T^      111 I      I   I  1 1  iiii|      I   I  I  I  iiiij      I   I  1 1  iiiij      I    I  I  I  iiiii — I   I  1 1  mil — I   I  M  iiiii — I    I  I  I  iiiii — I   I  I  mill — r 


iiil — I   I  I  Mini — ml — I   I  I  mill — il I   I  I  iiiiil I   I  I  mill     I  I  I  iiiiil     I  I  I  iiiiil     I  I  I I     I  I  r 

1012         io»3         10'*         1016         1016         iQn         iQie         IQi^         lO^o         lO^i 

M/Mg 


Figure  5-3.  The  bias  parameter  plotted  as  a  function  of  halo  mass.  Most 
halos  fall  into  the  mass  range  M/Mq  =  W^  -  W^,  where  the 
bias  is  near  unity. 
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Figure  5-4.  Dimensionless  nonlinear  power  spectra  using  the  NFW  halo  pro- 
file. The  sohd  hne  uses  the  Press-Schechter  mass  function,  the 
dashed  line  the  Sheth-Tormen  mass  function.  The  dotted  line 
shows  the  linear  power  spectrum  for  comparison. 
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Figtire  5-5.  Dimensionless  nonlinear  power  spectra  using  the  Moore  halo  pro- 
file. The  solid  line  uses  the  Press-Schechter  mass  fimction,  the 
dashed  line  the  Sheth-Tormen  mass  function.  The  dotted  line 
shows  the  hnear  power  spectrum  for  comparison. 
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Figure  5-6.  Dimensionless  nonlinear  power  spectra  using  the  NFW  profile  for 
large  mass  halos  and  the  Moore  profile  for  smaller  mass  halos. 
The  solid  line  uses  the  Press-Schechter  mass  function,  the  dashed 
line  uses  the  Sheth-Tormen  function.  The  dotted  line  shows  the 
linear  power  spectrum  for  comparison. 
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Figure  5-7.  Dimensionless  nonlinear  power  spectrum  compared  to  linear  ver- 
sion. The  solid  line  is  the  nonlinear  spectrum,  using  the  mixed 
halo  profile  and  the  ST  mass  function.  The  dotted  line  repre- 
sents the  linear  power  spectrum.  The  short  dashed  line  is  the 
one-halo  contribution  to  the  nonlinear  spectrum  and  the  long 
dashed  hne  is  the  contribution  from  the  halo-halo  term. 
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where  q  —  v{l  +  S)  is  the  density  weighted  pecuhar  velocity,  which  has  the  Fourier 
transform 

q  =  v  +  j  -^^v{k')5{k-k').  (5.24) 

The  angular  two-point  correlation  function  of  the  CMB  temperature  fluctua- 
tions is  defined  by  equation  [4.2].  The  power  spectrum,  Cg,  can  be  calculated  as  a 
projection  of  the  three-dimensional  power  spectrum.  This  result  is  due  to  Kaiser,^^ 
and  is  simply  the  Fourier  version  of  a  projection  originally  due  to  Limber.^'^ 

For  the  density  weighted  velocity,  Ma  and  Fry  arrive  at  an  expression  for  the 
angular  power  spectrum  for  the  kSZ  effect, 


>   ■■  1 


where  P^^  is  the  power  spectrum  of  q  along  the  line  of  sight  vector.  The  integration 
variable,  x,  is  simply  the  co-moving  distance  to  an  object  of  redshift  z.  In  a  model 
with  matter  and  cosmological  constant,  but  no  curvature,  the  relationship  between 

X  and  z  is  given  by         '  ■     ' 

dz 
dx  = ,  ^    ■  (5.26) 

Figure  [5-8]  shows  this  distance  plotted  against  {l-\-  z). 

One  very  interesting  aspect  of  a  projection  such  as  this  is  that  modes  parallel  to 
the  wave- vector  k  do  not  contribute,  as  the  crests  and  troughs  cancel  each  other  when 
collapsed  on  top  of  one  other."*^  (In  actuality  this  cancellation  is  only  approximate 
for  long  wavelengths,  or  large  angular  scales,  but  is  very  near  total  at  the  small 
angular  scales  we  are  interested  in.)  This  imphes  that  what  we  really  need  to 
find  is  Pqx  ~  2Pg-y.  The  factor  of  two  comes  from  the  fact  that  there  are  two 
independent  directions  perpendicular  to  fc.  It  also  implies  that  there  is  no  linear 
perturbative  contribution  to  this  eff'ect,  as  in  linear  theory  (without  the  rotational 
mode)  q  =  v{k)  oc  k.  It  is  the  second  order  contribution,  called  the  Ostiker-Vishniac 
(OV)  effect,  that  we  will  look  at  in  the  next  section. 
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Figure  5-8.  The  relationship  between  redshift  and  distance  to  an  object  in  a 
zero-curvature  universe. 
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The  result  obtained  by  Ma  and  Fry  for  Pq±  is 


9± 


*^)  =  /s 


(l-/i'2)P^,(|fc-fc'|)P,,(A;') 


^^      ^"]^Ps.{\k-k'\)Ps.ik') 


I 


\k-k' 
d^k'  d^k" 


v/l-//'Vl-/^"^cos((/)'  -  0") 


(27r)3  (27r)3 

where  fj/  =  k  ■  k  ,  fj,"  =  k  •  k" ,  and  the  power  spectra  are  defined  by 

{S{ka)S{kb))    =    {2nf5D{ka  +  h)Pssika)  (5.28) 

{V'{ka)v^{h))      =      {2TTfSDika  +  h)UXPvv{ka)  (5.29) 

{S{ka)v'{h))    =    {2-Kf5n{ka  +  kb)k\Ps,{k,)  (5.30) 

{6{ka)~5{k,)v\K)v^{kd))c     =      {2'Kf5D{Y.k^)k%P65vv{kaM.K,kd).    (5.31) 

The  subscript  c  in  the  last  equation  refers  to  the  connected,  or  irreducible,  moment. 
It  is  important  to  note  that  equations  [5.29],  [5.30],  and  [5.31]  all  depend  on  the 
assumption  of  v  oc  fc,  a  condition  we  will  later  modify. 

5.5     Ostriker-Vishniac  Angular  Power  Spectrum 

Calucluation  of  the  angular  power  of  CMB  fluctuations  in  the  quasi-linear 
regime,  using  second-order  perturbation  theory,  produces  the  Ostiker-Vishniac  spec- 
trum. 

In  this  regime,  we  use  first  order  peculiar  velocities  and  density  perturbations 
to  get  the  second  order  power  spectra.  Taking  the  Fom-ier  transform  of  equation 
[2.22],  we  have 

v{k)  =  ^~5k  (5.32) 

A/ 
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so  that 

n™    =     [%)'pS'  (5.33) 

P^    =     {^i)p&'-  (5.34) 

Here  the  superscripts  refer  to  the  hnear  order  power  spectra.    Putting  these  in 
equation  [5.27]  we  have 


p(2) 


Using  the  hnear  power  spectrum  given  in  section  [5.2]  and  inserting  this  in  equation 
[5.25]  produces  the  spectrum  plotted  in  figure  [5-9].  In  that  figure  we  have  integrated 
equation  [5.25]  from  the  present,  2;  =  0,  to  several  different  possible  redshifts  of 
reionization,  Zr-  The  OV  spectrum,  and  the  more  refined  nonlinear  spectrum  of  the 
next  section,  are  sensitive  to  Zr  and  thus  observations  of  this  effect  will  put  good 
constraints  on  the  ionization  history  of  the  universe. ^^ 

5.6     Nonlinear  Angular  Power  Spectrum 
In  the  nonhnear  regime  Ma  and  Fry  suggest  an  approximate  form  for  P^x  which 
is  vahd  for  high  k,  the  region  of  interest  in  a  nonhnear  analysis.  They  use 

P,±{k)  =  IJ  0j^Pss{k)PUk')  =  l{v')Pss{k).  (5.36) 

In  figure  [5-10]  we  plot  the  angular  spectrum  obtained  from  this  form,  using  the 
one-halo  power  spectrum  given  in  section  [5.3].  We  use  a  value  of  (u^)  ^^  ^  200km/s 
for  the  velocity  dispersion.  The  key  feature  is  that  the  OV  spectrum  underestimates 
the  power  at  high  ts.  In  figure  [5-11]  we  show  the  nonlinear  kSZ  spectrum  as  it 
compares  to  the  expected  spectrum  from  the  thermal  SZ  effect.  Note  that  for  very 
high  i  the  kSZ  effect  is  dominant. 

One  factor  that  is  not  taken  into  account  in  this  model  is  the  bias  on  small 
scales  between  baryonic  matter  and  dark  matter.  The  SZ  effect  is  due  to  scattering 
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Figure  5-9.  The  second-order,  perturbative  Ostriker-Vishniac  angular  power 
spectrum.  The  soUd  hne  is  the  OV  spectrum,  calculated  using 
a  redshift  of  reionization  of  20,  the  dashed  line  is  the  same  with 
reionization  redshift  of  6,  and,  for  reference,  the  dotted  line  is 
the  primary  anisotropy  spectrum. 
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off  free  electrons  which  may  (in  fact,  most  certainly  do)  cluster  differently  than  the 
surrounding  dark  matter  halo.  Ma  and  Fry  do  address  this  by  also  plotting  the 
spectrum  using  a  somewhat  different  profile  function  based  on  the  modeling  of  hot 
gas  in  galaxies  and  clusters.  This  generally  causes  the  power  to  be  lower,  but  still 
greater  than  in  the  linear  model. 

Before  moving  on  to  the  next  chapter  and  the  effect  of  rotational  velocities  on 
these  power  spectra,  we  note  that  we  have  not  considered  the  possibility  of  non- 
imiform  reionization  in  these  calculations.'*^- ^°  Santos  et  al.  calculate  that  such 
"patchy"  reionization  would  lead  to  significantly  higher  values  for  the  temperature 
power  spectrum,  and  in  fact  would  likely  swamp  the  effect  of  density  perturbations. 
Since  httle  is  known  about  the  actual  inhomogeneity  of  reionization,  we  have  chosen 
to  assume  it  uniform  for  all  that  follows. 
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Figiire  5-10.  The  nonlinear  angular  power  spectrum,  based  on  the  model  of 
Ma  and  Fry.  (solid  line)  The  dashed  Une  is  the  OV  spectrum 
and  the  dotted  hne  is  the  primary  anisotropy  spectrum.  Both 
the  nonhnear  and  OV  spectra  are  evaluated  to  z^  =  20. 
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Figure  5-11.  The  solid  line  represents  the  nonlinear  kSZ  effect  spectrum, 
while  the  dotted  line  represents  the  thermal  SZ  spectrum.  (Plot 
for  this  spectrum  was  obtained  from  C.-P.  Ma.)  The  dashed  line 
is  the  primary  spectrum. 


CHAPTER  6 
SZ  POWER  SPECTRUM  INCLUDING  ROTATIONAL  MODES 

6.1     Introduction 

In  this  chapter  we  calculate  the  angular  power  spectrum  of  temperature  an- 
isotropies  in  the  CMB  due  to  the  kinetic  SZ  effect  with  the  inclusion  of  rotational 
velocity  modes.  Before  beginning  our  discussion,  we  note  that  attempts  have  been 
made  to  include  rotation  in  a  calculation  of  the  kinetic  SZ  effect.  Cooray  and  Chen 
(2002)^^  calculate  the  angular  power  spectrum  of  the  rotational  kSZ  effect  using  a 
perscription  for  the  halo  rotations  based  on  work  by  Bullock  and  collaborators.^^ 
Importantly,  however,  this  is  angular  momentum  gained  through  events  such  as 
halo  interactions  and  is  different  from  the  rotational  modes  that  would  be  present 
due  to  conservation  of  initial  angular  momentum.  Cooray  and  Chen  find  a  power 
spectrum  contribution  several  orders  of  magnitude  smaller  than  the  contributions 
seen  in  the  previous  chapter,  and  in  fact  they  anticipate  this  result  given  the  accuracy 
with  which  Ma  and  Pry's  spectrum  matches  the  N-body  data.  Of  course  these  N- 
body  simulations  do  not  include  the  initial  rotational  modes,  so  they  would  not  be 
expected  to  include  the  contribution  we  calculate  here. 

6.2     First-Order  Results 

We  now  include  the  effect  of  rotational  velocities  in  the  SZ  angular  spectrum  by 
following  the  procedure,  due  to  Ma  and  Fry,  laid  out  in  the  last  chapter.  Remem- 
bering that  we  are  looking  for  the  power  spectrum  of  the  density  weighted  pecuhar 
velocity,  we  had  written  the  Fourier  transform  of  this  quantity  as 


q{k)  =  v{k)  +  J  ^,v{k')~5{k  -  k').  (6.1) 
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We  have  previously  ignored  the  contribution  from  the  first  term  above,  because  for 
the  growing  mode  v  <x  k.  The  first  order  rotational  mode,  however,  is  perpendicular 
to  k  so  we  have  the  simple  relationship,  Pg±  =  {ai/a)P^^^^,  where  P^^^^  is  the 
power  spectrum  of  the  rotational  velocity  at  decoupling  and  ai  is  the  scale  factor  at 
decoupling. 

We  expect  the  rotational  velocity  dispersion  at  decoupling  to  be  approximately 
equal  to  the  longitudinal  dispersion,  which  is  approximately  20  km/s.  This  means 


I 


(27r) 


P^.Sk')  =  {(v^r)  ^  400(km/s)2.  (6.2) 


We  can  determine  the  fimctional  dependence  of  the  initial  spectrum  on  the 
wavenumber  by  considering  the  mass  contained  within  the  horizon  at  a  time  t. 

Mh  =  Porl  =  Poictnf  =  (67rGpo4)(^)-  (6.3) 

We  write  it  in  this  manner  because  the  solution  to  equation  [2.9],  for  zero  curvature, 

is  Q-kGpqI^  =  1.  This  means  that  tn  oc  Mh  oa  k~^.  Finally,  we  note  that 

Aik)  =  eP{k)  oc  (^)'  a  (^-f\       oc  {k-'Y'^  =  k-\  (6.4) 

which  means  that  Py^yj_{k)  =  Ak'"^ .  We  see  that  we  need  a  low  k  cutoff  for  the 
integral  above  to  converge.  We  choose  a  cutoff  of  A;  =  10~^(/i/Mpc).  This  gives  us 
the  value  ^  =  3  x  10"^^/i(Mpc)"^(km/s)2.  Inserting  this  spectrum  into  equation 
[5.25]  gives  us  the  angular  spectrum  of  CMB  temperature  fiuctuations.  The  result 
is  given  in  figure  [6-1].  As  expected,  the  contribution  is  many  orders  of  magnitude 
below  other  contributions  because,  to  linear  order,  the  rotational  mode  decays  with 
the  expansion  of  the  universe. 

6.3     Non-Linear  Results 
Seeing  that  the  first  order  rotational  contribution  is  negligible,  we  look  to  the 
second  term  in  equation  [6.1]. 
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Figure  6-1.  Angular  power  spectrum  of  temperature  variation  on  the  CMB 
due  to  the  first  order  rotational  velocities. 
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Before  we  do  so  we  need  a  way  to  describe  the  rotational  velocity  field.  We  find 
this  by  imposing  angular  momentum  conservation  on  a  collapsing  object.  As  we 
saw  previously,  collapsed  objects  have  density  profiles  as  given  in  equations  [5.11] 
and  [5.12].  We  consider  the  mass  contained  within  an  initial  co-moving  radius  x, 
in  a  region  with  a  nearly  homogeneous  density.  (The  "nearly"  indicating  the  slight 
overdensity  which  will  lead  the  region  to  collapse.)  We  denote  the  final  co- moving 
raxiial  coordinate  y  and  write  the  mass  as 


M  =  — pciflfa;^  =  poa^ 


y 


rv 
y^+  /    dy'4Try'\u{y/ys) 
Jo 


(6.5) 


which  is  just  expressing  the  fact  that  we  are  superimposing  the  halos  onto  a  imiform 
backgrotmd.  This  will  allow  us  to  define  the  mapping  from  xtoy  for  a  given  u{y/ys). 
For  the  Moore  profile  {u{x)  =  l/[x^-^(l  +  x^-^)])  we  find 


—po,iaiX   =  poa 


^y'  +  y!^sjn[i  +  {y/y,f^] 


(6.6) 


but  since  po  oz  a  ^,  po.iof  =  poo.^  and  we  have 


x'  =  y'  +  2S^y'jn[l  +  {y/ysf/'] 


(6.7) 


As  our  interest  is  in  the  small  y  regime,  we  find  the  asymptotic  behavior  at  low  y 
to  be 

x'  =  25a{yysf'^.  (6.8) 

Figiire  [6-2]  shows  this  behavior. 

For  the  NFW  profile  {u{x)  =  \/[x{\  +  xf])  we  find 

'^-{y/ys) 


x^  =  y'  +  3Sayi 


A  +  iy/ys) 


+  ]n[l +  {y/y,)] 


which  asymptotes  to 


x^  =  35ay 


(6.9) 


(6.10) 


for  low  y.  (See  figure  [6-3].) 
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Figure  6-2.  Mapping  of  the  initial  co-moving  radius,  x,  to  final  co-moving 
radius,  y  for  a  Moore  profile  cluster.  The  values  5a  =  1/2  and 
Vs  =  2  have  been  used.  The  dotted  line  represents  the  limiting 
behavior  at  low  y. 
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Figure  6-3.  Mapping  of  the  initial  co-moving  radius,  x,  to  final  co-moving 
radius,  y  for  a  NFW  profile  cluster.  The  values  6a  —  1/2  and 
t/s  =  2  have  been  used.  The  dotted  line  represents  the  limiting 
behavior  at  low  y. 
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We  can  express  conservation  of  angular  momentum  as 


ay        a 


(6.11) 


where 


A 


XM{y)    =    {l  +  2Sa{ys/yr\n[l  +  {y/y,f/']f\ 

'i-{y/ys) 


NFW 


(y) 


l  +  3Sa{ys/yy 


+  ln[l  +  (y/ys)] 


1/3 


(6.12) 
(6.13) 


l^  +  iy/ys) 

Now,  returning  to  the  definition,  q  —  {I  +  6)v,  we  see  that  in  the  nonhnear 
regime,  where  the  second  term  above  will  be  dominant,  we  have,  for  the  rotational 
velocity  mode,  q  =  Svt  =  {ai/a)SXvo^f  We  then  make  the  definition  5'  =  {ai/a)SX 
and  write  q  =  S'vo^f  This  may,  at  first,  seem  like  an  odd  thing  to  do,  but  we  will 
find  that  we  just  absorb  the  collapse  factor  into  the  halo  density  profile  to  produce 
a  new  profile,  making  the  calculation  very  similar  to  that  in  the  last  chapter. 

The  second  moment  of  the  components  of  q  obey 

{q\k,)q^ik2))    =    {2nfSo{kx  +  k2)\{v\vi) 


+ 


/J6U      JAUI  ^  . 

^^(*^(fc')^(fci  -  k'W{k")5{k,  -  k"))] 


=    (27r)%(fci  +  fc2) 


{vW,) 


(M^  (2^  [(^^(fc')'^(fci  -  k')){v^{k")5{k,  -  k")) 
+{v\k')S(ki  -  k")){v\k")S{k2  -  k')) 


+{v'ik')v^ik")){6{k,  -  k')5{k2  -  k"))  +  {), 


(6.14) 


In  the  nonhnear  regime  behavior  is  dominated  by  the  terms  of  higher  order,  so  we 
will  neglect  the  first  term  above,  which  is  the  lowest  order  term.  The  second  term 
above  is  zero  by  our  assumption  of  isotropy.  The  last  term  is  the  irreducible  fourth 
order  term.  If  we  had  kept  the  collapse  factor  with  the  rotational  velocity,  we  would 
expect  this  irreducible  term  to  be  the  term  of  interest,  as  then  the  velocity  would 
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couple  strongly  to  the  density.  However,  as  we  have  associated  the  collapse  with  the 
density,  we  now  expect  the  important  term  to  be  similar  to  the  important  term  for 
the  longitudinal  mode.  Thus  the  irreducible  fouth  moment  will  vanish,  just  as  in 
Ma  and  Fry,  and  we  drop  it. 

We  now  write  the  transformed  velocity  as 

v{k)  =  vii{k)  +  v±{k)  =  kvii{k)  +  iv±{k),  '  (6.15) 

where  i  represents  a  unit  vector  perpendicular  to  k.  This  term  could  be  written, 

e  =  ^^-  -  ^^^'^^^  (6.16) 

v± 
Using  this  velocity  we  have  the  general  definitions, 

{S{k,)~S{k2))    =    {2TrfSD{k,  +  k2)Pss{h)  (6.17) 

{v\k,)v^{k2))    =    ([t)i|(fei)  +  t;i(fci)][t)j(fe2)  +  t>i(fc2)]) 

=      {V\{k,)vl{k,))  +  {V\{k,)v^^ik2)) 

=    {2TTfSo{kr  +  k2)[k^MPvn{ki)  +  i\iiP,^vAki)]   (6.18) 
{S{kr)v'{k2))    =    {S{k,)v\{k2))  +  {S{k,)vi{k2)) 

=    {2n)%{ki  +  k2)[k^2Psv,{ki)+iiPs,Ak,)],  (6.19) 

The  cross  terms  in  the  second  expression  vanish  because  the  components  are  inde- 
pendent. The  last  term  in  equation  [6.19]  is  zero  by  isotropy. 

Now,  remembering  that  wherever  we  have  a  S  associated  with  a  rotational 
velocity,  we  will  replace  it  with  a  S',  and  we  will  replace  the  velocity  with  Vo,t,  we 
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have 

P-(k)    =     f^f^ 
'^''^   '  J    (27r)3(27r)3 

X  [{27rfSo{k  -k'-  fc")[fe"fe"P.„.|,(A:')A5(fc")  +  e"e"'P^^,^{k')PMk") 
+k"k"'Ps,^^  (A;')F<5.||  (A;")  +  fc^f'^P^^,,  {k')Ps',^  {k") 
+i"k"'Ps>.,^,{k')Ps,^^{k")  +  £'^e"^P,,,^_,(A;')P5'.,,o(A;")]].  (6.20) 

Note  that  e  ,  e  refer  to  unit  vectors  perpendicular  to  k  and  k  ,  respectively.  We 
now  calculate  the  quantity  of  interest  using 

P<,±(fc)  =  2f  7^PJJ(A:).  (6.21) 

The  inner  products  of  interest  are 

f^^fc^fe''    =    {l-fi'^y^  (6.22) 

f7^fcV    =    {fi'{l-fiy/^)u'^  (6.23) 

ri'k'k"'  =  (l-/i'2)i/2^'^.fe"  (6.25) 

fyfe'e"^  =  (l-/i'2)i/2^.£"  (6.26) 

f^f'^fe"^  =  //V'7-fc"  (6.27) 

f^f'^e"^  =  /iV'7-e",  (6.28) 

where  /i'  is  defined  as  before  and  v'  =  cos(0').  In  equations  [6.25]  and  [6.27]  we  have 
the  expression,  7  •  fc  ,  which  we  evaluate  by  noting  that  in  the  next  step  we  will 
integrate  over  the  delta  function  in  equation  [6.20],  and  thus  replace  k  with  k  —  k. 
Since  fc    ^  (fc  —  fc  )/|fc  —  fc  |  and  7  •  fc  =  0,  we  end  up  with 

^.yfc"r    =    ^i|^^  (6.29) 


.   ♦    'I 
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Additionally,  in  equations  [6.26]  and  [6.28],  7  •  i"  becomes  7  •  i'"  where  e'"  is  a  unit 
vector  perpendicular  to  fc  —  fc'. 
Putting  these  in,  we  have 

+(1  -  t^'Y'u'y  ■  e"'Ps.,ik')Ps>.,,{\k  -  fc'l) 

u'(l  -  u'^V/'^iy''^k' 
^\k^\ ^^'-x.o(^')A.„(|fc  -  fc'l) 

+/xV'7  •  f"'P,,,,  „(A:')Py.,,o(|fc  -  fc'l)].  (6.31) 

Now,  noting,  as  per  Ma  and  Fry,  that  we  are  interested  in  high  k,  beyond  the  k'  peak 
of  the  above  integrals,  we  drop  terms  of  order  k'/k  and  perform  the  integrations  over 
/i'  and  u'.  First,  we  notice  that 

7  •  i'"  ^  7  .  f  =  1,  (6.32) 

in  this  hmit.  To  perform  the  angular  integrations,  we  note  that 

/  (fk'  cos2(0')/(/i',  k,)    =  n  f   dn'  r  dk'k'^fie',  k') 

=  \ld'k'f{9',k'),  (6.33) 

fd^k'cos{(f>')f{^i',k')    =  0,  (6.34) 

Jd'k'{l-f,'')f{k')    =  2n^Jdk'k"f{k') 

=  IJd'k'fik'),  (6.35) 

Jd^k'fi'^fik')    =  2tt^  fdk'k'^fik') 

=  ^Jd'k'fik'),  (6.36) 

Jd'k'fi'{l-f,y/'fik')    =  0.  (6.37) 
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Putting  this  together  gives  us  our  result, 


/d^k'    [2  1 

j^s  hPnn{k')Pss{k)  +  -Pl,,{k')PMk) 


(6.38) 


The  first  term  is,  of  course,  the  term  introduced  last  chapter,  due  to  Ma  and 
Fry.  It  is  the  second  term  we  calculate  here. 

First,  we  must  calculate  Ps's'{k).  As  we  hinted  at  earher,  the  collapse  factor, 
A,  simply  multiplies  the  halo  density  profile,  giving  us  a  new  profile,  w{y/ys)  = 
Xu{y/ys).  Since  we  are  ultimately  interested  in  behaviors  at  low  y,  or  equivalently, 
high  k,  we  will  use  the  low  y  asymptotic  form  of  A.  For  the  Moore  profile, 

It  is  the  Fourier  transform  of  this  function  that  we  will  need.  Figure  [6-4]  shows 
this  transform  as  compared  to  the  transform  of  the  full  w.  The  discrepancy  is  is  the 
relatively  unimportant  small  k  region. 

Figure  [6-5]  shows  the  transform  of  w  and  an  analytical  fit,  which  we  will  use 
to  increase  computational  efficiency.  This  fit  is  given  by 

For  the  NFW  profile,  we  have 

This  transform  and  the  transform  of  using  the  full  version  of  A  are  shown  in  figure 
[6-6]. 

The  analytical  fit  of  the  transform  of  this  profile  is 

Figure  [6-7]  shows  this  fit  compared  to  the  transform. 
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0.001        0.01  0.1 


Figure  6-4.  The  Fourier  transform  of  the  Moore  w{y).  The  sohd  hne  rep>- 
resents  the  use  of  the  full  version  of  A,  while  the  dotted  line 
represents  the  use  of  the  high  k  asymptotic  form  of  A.  In  both 
versions  we  set  Sa  =  1/2. 
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Figure  6-5.  The  Fourier  transform  of  w{y)  for  the  Moore  profile  using  the 
asymptotic  form  of  A.  The  solid  line  represents  the  actual  trans- 
form and  the  dotted  line  represents  the  fit  given  in  equation 
[6.40] 
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Figure  6-6.  The  Fourier  transform  of  the  NFW  w{y).  The  solid  Hue  rei>- 
resents  the  use  of  the  full  version  of  A,  while  the  dotted  line 
represents  the  use  of  the  high  k  asymptotic  form  of  A.  In  both 
versions  we  set  da  =  1/3. 
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Figure  6-7.  The  Fourier  transform  of  w(y)  for  the  NFW  profile  using  the 
asymptotic  form  of  A.  The  soUd  hne  represents  the  actual  trans- 
form and  the  dotted  line  represents  the  fit  given  in  equation 
[6.42] 
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We  now  simply  calculate  Ps's'ik)  just  as  we  did  before,  in  equation  [5.19],  with 
the  substitution  of  w  for  u.  Figure  [6-8]  shows  Ps'S'ik),  with  the  factor  of  {ai/af 
left  out,  as  compared  to  Ps5{k).  Both  the  overall  magnitude  and  the  shape  are,  of 
coiurse,  markedly  different. 

Figure  [6-9]  shows  the  difference  in  Ps's'{k)  when  the  different  halo  profiles  are 
used.  Notice  that  use  of  the  NFW  profile,  rather  than  the  Moore  profile,  leads  to  a 
difference  of  about  a  factor  of  ten  in  the  spectrum  for  high  k. 

We  now  have 

P,^  =  \{vlt)P5'6'{k).  (6.43) 

We  will  use  a  value  of  {vQ^Y/'^  =  20km/s,  as  we  did  in  the  first  section  of  this  chapter. 

Putting  this  into  equation  [5.25],  we  find  the  result  given  in  figures  [6-10]  and 
[6-11]. 

Figure  [6-12]  shows  how  the  nonlinear  rotational  velocity  contribution  compares 
to  the  linear  contibution.  This  shows,  in  stark  fashion,  that  this  study  was  well 
motivated.  i 


...  -t. 
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Figure  6-8.  Plots  of  Fj'5'(A;),  the  solid  line,  and  P(55(fc),  the  dashed  line. 
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Figure  6-9.  P^-^'  {k)  plotted  using  the  NFW,  dashed  Une,  Moore,  dotted  hne, 
and  mixed,  soUd  hne,  halo  profiles. 
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Figure  6-10.  Angular  power  spectrum  of  temperature  fluctuations  in  the 
CMB  due  to  nonlinear  rotational  velocity  contribution. 
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Figure  6-11.  The  solid  line  shows  the  angular  spectrum  due  to  the  rotational 
velocity  mode.  The  nonlinear  spectrum  due  to  the  longitudi- 
nal mode,  dotted  line,  and  the  primary  anisotropy  spectrum, 
dashed  line,  are  shown  for  comparison. 
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Figure  6-12.  The  angular  power  spectrum  of  CMB  temperature  fluctuations 
due  rotational  velocity  modes.  The  top  hne  is  the  nonlinear 
contribution,  the  bottom  hne  the  linear  contribution. 


CHAPTER  7 
CONCLUSION 

Our  goal  in  this  dissertation  has  been  to  include  the  effects  of  an  initial  ro- 
tational velocity  mode  in  large-scale  structure  calculations.  We  have  retained  this 
mode  in  perturbation  theory  and  we  found  that,  although  it  decays  with  time  in 
linear  theory,  at  higher  orders  it  couples  to  the  longitudinal  mode  and  plays  a  non- 
trivial  role  in  quasi-linear  theory. 

We  expected  that  the  rotational  velocity  mode  might  play  an  even  more  signifi- 
cant role  in  the  nonhnear  regime.  Specifically,  we  beheved  that  the  inclusion  of  this 
mode  in  a  calculation  of  the  kinetic  Sunyaev-Zeldovich  effect  might  prove  significant. 
The  results  of  last  chapter  indicate,  however,  that  the  inclusion  of  such  modes  is  not 
significant  on  any  angular  scales  available  to  observation  now  or  in  the  near  future. 

This  fact  provides  evidence  for  the  vahdity  of  the  assumption  typically  made 
in  previous  calculations  that  the  rotational  mode  is  largely  unimportant  in  most 
cosmological  situations.  It  is  important,  however,  to  note  that  the  inclusion  of 
rotational  modes  in  the  nonlinear  calculation  leads  to  temperature  fluctuations  many 
orders  of  magnitude  higher  than  those  in  the  linear  regime.  It  seems  that  the  success 
of  the  a  priori  exclusion  of  initial  rotational  modes  was  fortuitous. 

It  is  important  to  note  that  there  are  many  possible  ways  to  improve  upon 
our  results.  One  important  factor  is  that  the  halo  profiles  used  to  build  the  halo 
model  are  the  result  of  N-body  simulations  which  do  not  include  the  initial  rotational 
mode.  It  is  certainly  conceivable  that  the  profile  in  the  central  regions  of  halos  would 
be  considerably  different  if  this  mode  were  included  in  the  simulations.  Analytical 
calculations  which  allow  for  the  inclusion  of  angular  momentum  do  indicate  that  a 
shallower  inner  profile  is  hkely.^^  One  would  expect  that  if  a  numerical  simulation 
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were  performed  where  the  initial  rotational  velocities  were  retained,  such  a  profile 
would  be  seen.  Having  said  this,  we  note  that,  as  suggested  by  the  fairly  small  profile 
diflFerences  seen  in  figure  [6-9]  even  a  somewhat  different  inner  profile  would  not  be 
likely  to  produce  an  observable  effect  on  the  angular  power  spectrum  of  temperature 
fluctuations. 

Besides  the  aforementioned  inclusion  of  rotational  modes  in  N-body  simulations, 
there  are  many  ways  in  which  the  work  begun  in  this  dissertation  could  continue. 
First,  the  role  of  the  rotational  modes  in  quasi-linear  theory  should  be  further  ex- 
plored. We  have  shown  here  that  it  plays  potentially  important  roles,  but  not  fully 
explored  what  those  roles  might  be.  Additionally,  the  inclusion  of  rotational  veloc- 
ities will  also  have  an  impact  on  the  thermal  Sunyaev-Zeldovich  effect.  An  initial 
rotational  velocity  dispersion  should  contribute  to  the  overall  velocity  dispersion, 
and  thus  effective  temperature,  of  collapsed  objects.  ,.,.. 
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APPENDIX  A 
CONSERVATION  OF  ANGULAR  MOMENTUM  IN  AN  EXPANDING 

UNIVERSE 

We  begin  by  defining  the  angular  momentum  about  an  arbitrary  origin  as 

L=      d^rpr  x  v,  (A.l) 

which  we  can  express  in  comoving  coordinates  as 

L=      (fxpoa^{l  +  5)ax  x  v.  (A.2) 

We  want  to  look  at  the  time  derivative  of  the  angular  momentum.  Noting  that  the 
quantity  PqO,^  is  constant  in  time  (remembering  that  po  on  a~^)  we  have 


dL 
Itt 


(fxx  X  —  [(1  +  6)av].  (A.3) 


Using  the  fluid  flow  equations  from  chapter  2,  equations  [2.12]  and  [2.15],  this  be- 
comes 

dL 


dt 


=    -poo    /  d-^x 


I  d^x  \x  X  {(1  +  8)[{v  •  V)v  +  V0]  +  v(V  •  [(1  +  5)v])}  (A.4) 
=    -pqo^  I d^x\xx{{[V  ■{l  +  5)v\v)  +  {\  +  5)V4>}\.  (A.5) 


The  cross  product  of  x  with  the  second  term  vanishes  because  (j)  leads  to  a 
central  force.  Writing  the  remainder  in  component  notation,  we  have 


dLi 
~dt 


=      d^xeijkXj'Vi[{l-S)vkVi],  (A.6) 


which  we  can  write  as 


dLi 


=  /fe«»lV,x,(l+^K.,  +  (l+.)».„,V,x,l.  (A.7) 
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The  first  term  is  a  surface  term  and  is  thus  zero,  and  the  second  term  is  just 


dLi 


=  /  (fxSijkSijVkViil  +  (5)  =  0,  (A.8) 


as  it  is  just  proportional  to  the  cross  product  of  v  with  itself.  Thus  dL/dt  —  0, 
and,  as  we  would  suspect,  the  total  angular  momentum  in  an  expanding  universe  is 
conserved. 


APPENDIX  B 
SPHERICAL  COLLAPSE  AND  THE  DERIVATION  OF  5c 

A  reasonable  value  for  the  critical  density  contrast  required  for  collapse  of  a 

halo  can  be  obtained  through  a  simple  analysis  of  spherical  collapse  model.    We 

begin  by  writing  the  Newtonian  equation  of  motion  for  a  mass  shell  surrounding  a 

mass,  M. 


(fr         GM 


This  can  be  solved  using  the  parametric  approach, 


(B.l) 


r  =  ^(l-cos^)  (B.2) 

t  =  B{e-sme).  (B.3) 

Simple  substitution  will  verify  this  solution  and  provide  the  relation 

A^  =  GMB^.  (B.4) 

Figure  [B-1]  shows  the  radius  and  time  plotted  as  fimctions  of  the  parameter 

e. 

Now,  expanding  these  solutions  for  small  9,  corresponding  to  small  initial  over- 
densities,  we  have 
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Figure  B-1.  Plots  of  radius  and  time  as  functions  of  the  parameter  6.  t/A  is 
the  solid  line,  r/B  is  dashed. 
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Feeding  the  leading  order  solution  for  t  into  the  expression  for  r  we  find 


r{0) 


A 
2 


B 


2/3 


20  \b) 


■-{mf/\GMyi^ 


^-h 


'&f 


2/3' 


(B.7) 
(B.8) 


Remembering  that  p  —  M/(47rr^/3)  (this  is  actually  the  mean  density  within  a 
sphere  of  radius  r)  and  po  =  l/Girt'^G,  we  find 


P-  Po 


&t\ 


2/3 


(B.9) 


Po         20  \BJ  ■    ■<:- 

Equation  [B.2]  tells  us  that  total  collapse  occurs  at  ^  =  2n,  where  t  =  2ttB.  Putting 
this  in  we  get  the  result,  6c  ~  1.686.  The  reader  may  object,  as  ^  =  27r  is  clearly 
not  <C  1,  however,  all  this  is  saying  is  that  collapse  takes  place  when  Unear  theory 
predicts  an  overdensity  of  1.686. 

As  a  point  of  further  interest,  we  can  calculate  the  exact  density  rim  for 
this  spherical  collapse.   To  start  we  consider  the  energy  of  a  spherical  shell  at  its 

turnaround, 

GM(x)  '    --. 

E= r4-^,  ■  :^^  (B.IO) 


2A     ' 

as,  at  turnaround,  all  energy  is  potential  and  r  =  2 A.  Here  x  refers  to  the  co- moving 
radius.  We  also  note  that  E  oc  1/x  because  we  are  perturbing  an  Einstein-De  Sitter 
universe,  where  £■  =  0,  by  an  initial  overdensity  a  co-moving  distance  x  from  the 
shell.  When  combined  with  the  fact  that  M{x)  oc  x^  we  see  that  A  (x  x"^  and 
B  oc  x^/^.  Using  these  scale  relations  we  can  write 


^(1-cos^)       f  x'^\  1-cos^ 


Tta 


2A 


and 


X 


ta 


9/2 


X 


ta 


B 


TV 


B 


ta 


e-sm9' 


(B.ll) 


(B.12) 
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so  that 

r        /       TT       \  '    1  —  cos  9 

In  each  of  the  previous  equations  the  subscript  "ta"  indicates  a  quantity  for  the  shell 
at  turnaround  at  the  chosen  time  t.  We  compute  the  mass  density  run  using 

,  .      dM    1  ,„      ^ 


and  using  the  chain  rule, 

dM        dMH.T  /dr\'^ 

(B.15) 


dr        dx  dO  \dd )      ' 

and  the  dependences  we  found  above.  The  result,  when  divided  by  the  mean  density, 

becomes 

p(r)  _  9(^-sin^)2  \^      9sin6'(^-sin6»)^~^ 


Po        2(l-cos6')3 


4- 


(B.16) 


2    (l-cos^)2 

We  can  also  write  down  the  mean  density  within  the  sphere,  to  compare  with  the 
linear  prediction, 

P  ^    ^  a^ri^  _  ^^Gt^  _  9(^-sing)^ 

J,  ~  l^^""^^   -  2~^  -  2  {I -cos  or  ^^-^^^ 

Figure  [B-2]  shows  the  density  contrast  as  a  function  of  the  scaled  radius,  for  the 
mean  density,  as  above,  and  for  the  Unear  prediction  of  equation  [B.9].  Also  shown 
the  the  density  run.  None  of  these  are  expected  to  be  accurate  at  small  r,  where 
shell  crossing  causes  virialization. 
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Figure  B-2.  Density  contrasts  plotted  as  a  function  of  the  scaled  radius.  The 
dashed  line  is  the  exact  solution  for  the  mean  density  contrast 
within  the  radius.  The  dotted  line  is  the  linear  model  prediction. 
The  solid  line  is  the  actual  density  run. 


APPENDIX  C 
PHYSICAL  CONSTANTS  AND  COSMOLOGICAL  PARAMETERS 


Table  C-1.  A  List  of  constants  and  parameters  used  in  calculations. 


Parameter 

Symbol 

Value 

Units 

Speed  of  Light 

c 

299792.458 

km/s 

Mass  of  Proton 

nip 

1.672621  X  10-27 

kg 

Gravitational  Constant 

G 

6.673  X  10-2° 

km^kg-is-2 

Thompson  Cross-Section 

(7t 

6.65245  X  10-^5 

km2 

Fractional  Matter  Density 

^m 

0.3 

- 

Fractional  A  Density 

Oa 

0.7 

- 

Fractional  Baryon  Density 

^b 

0.05 

- 

Primordial  Helium  Fraction 

Y 

0.74 

- 

Ionization  Fraction 

Xe 

1.0 

- 

Initial  Power  Index 

n 

1.0 

- 

Dimensionless  Hubble  Constant 

h 

0.71 

- 
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APPENDIX  D 
CODE  FOR  CALCULATING  ANGULAR  SZ  SPECTRUM 

The  integration  elements  of  this  code  were  adapted  from  routines  given  in  Press 

et  al.^ 

/*  rvcl.c 

Calculates  the  rotational  contribution  to  Cl's 


-*/ 


#include  <math.h> 

#include  <stdio.h> 

#include  " ./cosmodef s.h" 

#define  FUNC(x)  (C+funk) (1 .0/(x))/((x)*(x))) 

#define  FUNCY(x)  ((*funk)(x)) 

#define  FUNT(x)  ((*funk)  (tan (x))* (1.0+ (tan (x))* (tan (x)))) 

#define  EPS  l.Oe-3 

#define  JMAX  14 

#define  JMAXP  15 

#define  K  5 

#define  NR.END  1 

#define  FREE_ARG  char* 

#define  ZMAX  20.0 

#define  MMAX  0 

#define  TMAX  30.0 

#define  LLMAX  6.0 

#define  ZZMAX  6.0 

#define  Do  0.31124 

#define  V2  400.0 

double  kay,  zltemp,  rad,  rsptemp,  kprimetemp; 
int  1; 

double  qromo(double  (*func) (double) ,  double  a,  double  b, 

double  (*choose) (double (*) (double) ,  double,  double,  int)) 
/♦Romberg  Integration  on  an  open  interval*/ 

{ 

void  polint(double  xa[] ,  double  ya[] ,  int  n,  double  x,  double  *y 

,  double  *dy) ; 

int  j; 

double  ss,   dss,   h[JMAXP+l] ,    s [JMAXP]; 

hCl]=1.0; 

for   (j=l;j<=JMAX;j++)    i 

s[j]=(*choose) (func,a,b, j) ; 
if   (j   >=  K)   { 

polint(&h[j-K] ,&s[j-K] ,K,0.0,&ss,&dss) ; 
if   (fabs(dss)    <=  EPS*fabs(ss))    return  ss; 
} 
h[j+l]=h[j]/9.0; 

} 

printfC'Too  many  steps  in  routine  qromoXn"); 
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exit(l); 

return  0.0:  ^    .,':.' 

double  midtan (double  (*fimk) (double) ,  double  aaa,  double  bbb,  int  n) 

double  x,tnin,sum,del,ddel,b,a; 
static  double  s; 
int  it , j  ; 

b=atan(bbb) ; 
a=atan(aaa); 
if    (n==l)   { 

return   (s=(b-a)*FUNT(0.5*(a+b)))  ; 
}  else  { 

for   (it=l,j=l;j<n-l;j++)    it  *=  3; 
tnin=it ; 

del= (b-a) / (3 . 0*tnm) ; 
ddel=del+del; 
x=a+0 . 5*del ; 
suin=0 . 0 ; 

for   (j=l;j<=it;j++)    { 
sum  +=  FUNT(x); 
X  +=  ddel; 
sum  +=  FUNT(x); 
X  +=  del; 

return   (s=(s+(b-a)*sum/tnm)/3.0) ; 

} 

double  midtanl (double  (*funk) (double) ,  double  aaa,  double  bbb,  int  n) 

double  x,tnm,sum,del,ddel,b,a; 
static  double  s; 
int  it,j; 

b=atan(bbb); 
a=atan(aaa) ; 
if  (n==l)  { 

return  (s=(b-a)*FUNT(0.5*(a+b))) ; 
}  else  { 

for  (it=l,j=l;j<n-l;j++)  it  *=  3; 
tnm=it ; 

del= (b-a) / (3 . 0*tnm) ; 
ddel=del+del ; 
x=a+0 . 5*del ; 
sum=0 . 0 ; 

for  (j=l;j<=it;j++)  { 
sum  +=  FUNT(x); 
X  +=  ddel; 
sum  +=  FUNT(x) ; 
X  +=  del ; 

return  (s=(s+(b-a)*sum/tnm)/3.0) ; 

} 

double  midtan4 (double  (*funk) (double) ,  double  aaa,  double  bbb,  int  n) 

double  x,tnm,sum,del,ddel,b,a; 
static  double  s; 
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int  it.j; 

b=atan(bbb) ; 
a=ataii(aaa) ; 
if   (n==l)   { 

return  (s=(b-a)*FUNT(0.5*(a+b))) ; 

for   (it=l,j=l;j<n-l;j++)    it  *=  3; 

tiiin=it ; 

del=(b-a)/(3.0*tnin); 

ddel=del+del ; 

x=a+0 . 5*del ; 

s\im=0 . 0 ; 

for   (j=l;j<=it;j++)    { 

smn  +=  FUNT(x); 

X  +=  ddel; 

sum  +=  FUNT(x); 

X  +=  del; 

return  (s=(s+(b-a)*sum/tnm)/3.0) ; 
} 
} 

double  midtanS (double   (*funk) (double) ,   double  aaa,   double  bbb,    int  n) 

{ 

double  x,tnin,sum,del,ddel,b,a; 
static  double  s; 
int  it , j  ; 

b=atan(bbb)  ; 
a=atan(aaa) ; 
if  (n==l)  { 

return  (s=(b-a)*FUNT(0.5*(a+b))) ; 

for  (it=l,j=l;j<n-l;j++)  it  *=  3; 

tnin=it ; 

del=(b-a)/(3.0*tnm); 

ddel=del+del ; 

x=a+0.5*del; 

sum=0 . 0 ; 

for  (j=l;j<=it;j++)  { 

sum  +=  FUNT(x); 

X  +=  ddel; 

sum  +=  FUNT(x); 

X  +=  del; 

return  (s=(s+(b-a)*sum/tnm)/3.0) ; 
} 
} 

double  midtan6 (double  (*funk) (double) ,  double  aaa,  double  bbb,  int  n) 

{ 

double  x,tnm,sum,del,ddel,b,a; 
static  double  s; 
int  it,j; 

b=atan(bbb) ; 
a=atcin(aaa) ; 
if  (n==l)  i 

return  (s=(b-a)*FUNT(0.5*(a+b))) ; 

for  (it=l,j=l;j<n-l;j++)  it  *=  3; 
tnm=it ; 


97 


del=  (b-a)  /  (3 .  0*tnin) ; 

ddel=del+del ; 

x=a+0 . 5*del ; 

simi=0 . 0 ; 

for  (j=l;j<=it;j++)  { 

sum  +=  FUNT(x) ; 

X  +=  ddel; 

sum  +=  FUNT(x) ; 

X  +=  del; 

return  (s=(s+(b-a)*sum/tnm)/3.0) ; 


double  midpnt (double  (*funk) (double) ,  double  a,  double  b,  int  n) 
/♦Integration  workhorse*/ 

{ 

double  x,tnm, sum, del, ddel; 
static  double  s; 
int  it ,  j  ; 

if  (n==l)  { 

return  (s=(b-a)*FUNCY(0.5*(a+b))) ; 
}  else  { 

for  (it=l,j=l;j<n-l;j++)  it  *=  3; 

tnm=it ; 

del= (b-a) / (3 . 0*tnm) ; 

ddel=del+del ; 

x=a+0.5*del; 

sum=0 . 0 ; 

for  (j=l;j<=it;j++)  { 

sum  +=  FUNCY(x) ; 

X  +=  ddel; 

sum  +=  FUNCY(x); 

X  +=  del; 

return  (s=(s+(b-a)*sum/tnm)/3.0) ; 


double  midpntKdouble  (*funk)  (double) ,  double  a,  double  b,  int  n) 
/♦Integration  workhorse*/ 

{ 

double  x.tnm, sum, del, ddel; 
static  double  s; 
int  it , j  ; 

if  (n==l)  { 

return  (s=(b-a)*FUNCY(0.5*(a+b))) ; 
}  else  { 

for  (it=l,j=l;j<n-l;j++)  it  *=  3; 

tnm=it ; 

del=(b-a)/(3.0*tnm) ; 

ddel=del+del ; 

x=a+0 . 5*del ; 

sum=0 . 0 ; 

for  (j=l;j<=it;j++)  { 

sum  +=  FUNCY(x); 

X  +=  ddel; 

sum  +=  FUNCY(x);  •  , , 

X  +=  del; 
} 
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return  (s=(s+(b-a)*sum/tnm)/3.0) ; 
} 
> 

double  midpnt2 (double  (*funk) (double) ,  double  a,  double  b,  int  n) 
/♦Integration  workhorse*/ 

{ 

double  x,tnin,sum,del,ddel; 
static  double  s; 
int  it, j; 

if  (n==l)  { 

return  (s=(b-a)*FUNCY(0.5*(a+b))) ; 

for  (it=l,j=l;j<n-l;j++)  it  *=  3; 

tnin=it; 

del= (b-a) / (3 . 0*tnm) ; 

ddel=del+del ; 

x=a+0.5*del; 

suin=0 . 0 ; 

for  (j=l;j<=it;j++)  { 

sum  +=  FUNCY(x) ; 

X  +=  ddel; 

sum  +=  FUNCY(x); 

X  +=  del; 

return  (s=(s+(b-a)*sum/tnm)/3.0) : 
> 
} 

double  midpnt4 (double  (♦funk) (double) ,  double  a,  double  b,  int  n)        S 
/♦Integration  workhorse*/ 

{  .: 

double  x.tnm, sum, del, ddel;  ■   , 

static  double  s;  '   ?;  v 

int  it.j;  v 

if  (n==l)  { 

return  (s=(b-a)+FUNCY(0.5^(a+b)));  -' 

for  (it=l,j=l;j<n-l;j++)  it  ♦=  3;  vV^:'.;' 

tnm=it;  "■'^  ■';•■?.■ 

del=(b-a)/(3.0^tnm);  ' 

ddel=del+del ; 

x=a+0.5*del; 

sum=0.0; 

for  (j=l;j<=it;j++)  { 

sum  +=  FUNCY(x); 

X  +=  ddel; 

sum  +=  FUNCY(x); 

X  +=  del;  :. 

return  (s=(s+(b-a)^sum/tnm)/3.0) ; 

} 

} 

double  *vector(long  nl,  long  nh) 

{  ■-   >■ 

double  ♦v;  .' 

v=(double  ♦)malloc((size_t)  ((nh-nl+l+NR_END)^sizeof  (double))) ;      ■'  '  .- 
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if  (!v)  { 
printf ("allocation  failure  in  vector\n"); 
exit(l);} 
return  v-nl+NR_END; 

void  free_vect or (double  *v,  long  nl,  long  nh) 

-C 

free((FREE_ARG)  (v+nl-NR_END)) ; 

void  polint (double  xa[],  double  ya[] ,  int  n,  double  x, 

double  *y,  double  *dy) 

/♦Polynomial  Extrapolation*/ 

{ 

double  *vector(long  nl,  long  nh) ; 

void  free_vector (double  *v,  long  nl,  long  nh.) ; 

int  i,m,ns=l; 

double  den,dif ,dift,ho,hp,w; 

double  *c,*d; 

dif =f abs(x-xa[l] ) ; 
c=vector(l,n) ; 
d=vector(l,n) ; 
for   (i=l;i<=n;i++)    { 

if   (   (dift=fabs(x-xa[i]))    <  dif)   { 

ns=i; 

dif=dift; 

c[i]=ya[i]  ; 

d[i]=ya[i]; 

} 
*y=ya[ns— ]  ; 

for   (m=l;m<n;m++)    { 
for   (i=l;i<=n-m;i++)    { 

ho=xa[i]-x; 

hp=xaCi+m]-x; 

w=c[i+l]-d[i]; 

if    (    (den=ho-hp)   ==0.0)   { 
printf ("Error  in  routine  polint\n") ; 
exit(l);} 

den=w/den ; 

d[i]=hp*den; 

c[i]=ho*den; 

*y  +=   (*dy=(2*ns  <   (n-m)   ?  c[ns+l]    :   d[ns~])); 

} 

free_vector(d, l,n) ; 

free_vector(c,l,n) ; 

double  ps  (double  x)   /*  This  is  the  k  dependent  peirt  of  the  linear 

power  spectrum  */ 

■c 

double  ans,  fac,  q,  p,  shape.  A,  delh,  omegamz,  omegalz,  f,  g,  D; 

delh=l . 94e-5* (pow(OMEGAM,  (-0 . 785-0 . 05* (log(QMEGAM) ) ) ) ) * 
(exp (-0 . 95* (INDEXN-1) -0 . 169* (INDEXN-1) * (INDEXN-1) ) ) ; 
A=delh*delh*pow (SOL/ (H* 100.0) ,  INDEXN+3)/(4.0*PI) ; 
shape=OMEGAM*H/ (exp(OMEGAB* (1 . 0+(l . 3/OMEGAM) ) ) ) ; 
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q=x/(shape*H) ; 
p=log(1.0+2.34*q)/(2.34*q); 

fac=1.0*(2.0*PI)*(2.0*PI)*(2.0*PI);   /*  To  account  for  difference 

between  Ma  and  Fry  */  - 

omegamz=OMEGAM/(OMEGAM  +  OMEGAL/(zltemp*zltemp*zltemp)) ; 

omegalz=1.0-omegamz;  ' 

f=pow(omegainz,  0.5714);  - 

g=omegamz/(f  -  omegalz  +  (1 .0+0.5*omegamz)*(l .0+0.014*omegalz)) ; 

/*  a  factor  of  2.5  left  out  of  g  and  of  Do*/ 

D=g/zltemp; 

ans= (D/Do) * (D/Do) *A*f ac*pow(x ,  INDEXN) *p*p/ (sqrt ( 1 . 0+3 . 89*q+16 . 1* 

16.  l*q*q+5. 46*5. 46*5. 46*q*q*q+pow(6.71*q,  4)));  ■*■_,:- 

return  ans:  ','?,'*?,; 


double  windowCdouble  x)  .. 

/*Fourier  transform  of  the  top-hat  window  function  */ 

{  ■   ■ 

double  win; 

win=3.0*((sin(x*rad))/(x*x*x*rad*rad*rad)-(cos(x*rad))  -  " 

/ (x*x*rad*rad) ) ; 

return  win;  '< 

double  grand  (double  x)  :  ■  •  ^^'^v 

return  (x*x*ps(x)*window(x)*window(x)/(2.0*PI*PI)) ;  '    :'  t 

}  ■•■■■::■.  ,-■,  :-■■_,:; 


double  cee (double  Inm) 
{ 

/*  concentration  factor  */ 
double  result; 

rad=pow( (3 . O*pow (10.0,  Inm) / (4 . 0*P1*AVDEN) ) 
,  (0.33333333333333)); 

result=qromo (grand,  0.0,  l.OeSO,  midtan6) ; 
if  (Inm  <  14.0)  { 

return  5 . 0*sqrt (result) ; 
}  else  { 

return  9. 0*sqrt (result) ; 

> 
} 


double  rs (double  Inm) 
i 

return  ((1.0/cee(lnm))*pow(3.0*pow(10.0,  Inm) 
/ (800 . 0*PI*AVDEN) ,   0 . 333333333333) ) ; 

double  stmf  (double  Inin) 


-^  •-•£\-iMr-i;^' 
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{ 

/*  Sheth-Tormen  mass  fxinction  */ 

double  dsidm,  deltaexp; 

double  result,  radtemp; 

double  invsigmaplus,  invsigmaminus ; 

double  resultplus,  resultminus ,  massfunc,  ma; 

ma=pow(10.0,    Inm) ; 

rad=pow((3.0*ma/(4.0*PI*AVDEN)) ,    (0.33333333333333)) ; 

result=qromo (grand,    0.0,    1.0e30,   midtaii4) ; 

deltaexp=ma* (1 . Oe-4) ; 

radtemp=rad; 

rad=pow( (3 . 0* (ma+deltaexp) / (4 . 0*PI*AVDEN) ) , 

(0.33333333333333)); 
resultplus=qromo (grand,  0.0,  1.0e30,  midtan4) ; 
invsigmaplus  =  1.0/(sqrt(resultplus)) ; 
rad=pow( (3 . 0* (ma-deltaexp) / (4 . 0*PI*AVDEN) ) , 

(0.33333333333333)); 
resultminus=qromo (grand ,  0.0,  1 . 0e30 ,  midtan4) ; 
invsigmaminus  =  1. 0/sqrt (resultminus) ; 
dsidm=(invsigmaplus-invsigmaminus)/(2.0*deltaexp) ; 
rad=radtemp ; 

massf unc=result*2 . 0*0 . 322*dsidm*AVDEN* (1 . 0+ 
1 . 0/(pow( (1 . 189*DC/sqrt (result) ) , 
0.6))) *sqrt ( 1 . 189* 1 . 189*DC*DC/ (result*2 . 0*PI) ) 
*exp (-DC*DC/ (0 . 707*2 . 0*result) )/ (ma) ; 
return  0.75*massfunc; 
} 

double  ftnfw (double  Inm) 
i 

/*  Algebraic  approximation  to  Fourier  transform 
of  NFW  halo  profile  */ 

double  q; 

q=kay*rs(lnm) ; 

return  (4.0*PI*(log(E+1.0/q)-log(log(E+1.0/q)/3.0)) 

/pow((1.0+pow(q,  1.1)), 

(2.0/1.1))); 

double  ftm  (double  Inm)  *■ 

{ 

/*  Algebraic  approximation  to  Fourier  transform  .  ,^ 

of  Moore  halo  profile  */  ,  .../ 

double  q;  >  v/   "  "m* 

q=kay*rs(lnm) ;  ,  .-  ■  ,   .    - 

return  (4.0*PI*(log(E+1.0/q)+0.25*log(log(E+1.0/q)))  -.' 

/(1.0+0.8*pow(q,  1.5)));  ;  '-r' 

}  "■  .  ■    /'-.-■-^CM 

double  ftmp(double  Inm)  '  ■■  ■'^; 


.:.:.- ;:^i'^ 


■^^ 
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/*  Algebraic  approximation  to  ft  of  Moore  profile  times 
collapse  factor  (high  k)  */ 

double  q; 

q=kay*rs(lnm) ; 

return  12.566*(log(E+1.0/pow(q,  0.3333))+0.25*log 
(log(E+1.0/pow(q.  0.33333) )))/pow(1.0+0.6*pow(q,  1.3),  0.77); 
} 

double  ftnfwp (double  Inm) 
{ 

/*  algebraic  approximation  to  ft  of  NFW  profile 
times  collapse  factor  (high  k)  */ 

double  q; 

q=kay*rs(lnm) ; 

return  12 . 566/ (1 . 0+0 . 6*q) ; 

} 

double  nfwd (double  Inm) 

/*  density  amplitude  for  NFW  profile  */ 

return  (200 . 0*cee (Inm) *cee (Inm) *cee (Inm) 

/  (3 .  0*  (logd  .  0+cee (Inm) )  -cee (Inm)  / (1 . 0+cee  (Inm) ) ) )  )  ; 

} 

double  md (double  Inm) 
{ 

/*  density  amplitude  for  Moore  profile  */ 

return  (100.0*cee(lnm)*cee(lnm)*cee(lnm) 
/logd. 0+pow(cee (Inm)  ,  1.5))); 

} 

double  integrand (double  zzl) 

i 

double  hofz; 

hofz=1.0/sqrt(0MEGAM*zzl*zzl*z2l  +  OMEGAL) ; 
return  hofz; 
} 

double  ex(double  zl)  /*  Calculating  comoving 

postion  as  a  function  of  z  */ 

{ 

double  ans2; 

if  (zKl. 00001) 
ans2=zl-l . 0-3 . 0*OMEGAM* (zl-1 . 0) * (zl-1 . 0) /4 . 0 ; 

ans2=(S0L/(H*100.0)) 
*qromo( integrand,  1.0,  zl,  midpnt4) ; 
return  ans2; 
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double  onehmintegrand (double  Inm) 
{ 

/*  Integrand  for  single  halo  power. 
Factor  of  mass  from  changing  to  dlnM*/ 

double  que,  rstempl,  mdtemp,  nfwdtemp,  ftmptemp, 
ftnfwptemp,  prefacm,  prefacnfw; 

rstempl=rs(lnm) ; 

que=kay*rstempl ; 

mdtemp=md(lnm) ; 

nfwdtemp=nfwd(lnm) ;      ■  .  '  ' 

ftmptemp=ftmp(lnm) ; 

ftnfwptemp=ftnfwp(lnm) ; 

prefacm=pow(2.0*mdtemp,  0.66666666) ; 

prefacnfw=pow(3.0*nfwdtemp,  0.66666666) ; 

if  (Inm  >  14.0)  { 

return  pref acnfw*pow(10.0,  lnm)*stmf (lnm)*pow(rstempl, 
6 . 0) *nf wdtemp*nf wdtemp*f tnf wptemp*f tnf wptemp ; 

}  else  { 

return  pref acm*pow( 10.0,  lnm)*stmf (Inm)* 
pow(rstempl ,  6.0) ♦mdtemp*mdtemp*f tmptemp*f tmptemp ; 

} 
} 

double  zintegrand (double  zl)  /*  Everything  to  be 

integrated  from  z=0  to  zr  */ 


i 


double  answer,  hover,  omegamz,  omegalz,  dxtodz; 

double  g,  D,  f,  pref,  extemp,  tauh,  tau,  expt,  lamafac; 

int  m,  t; 

extemp=ex(zl) ; 

kay=l/extemp ; 

zltemp=zl; 

answer=(1.0/3.0)*V2* 

qromo (onehmintegrand,  5.0,  18.0,  midpntl) ; 

dxtodz=S0L/(H*100.0*sqrt(0MEGAM*zl*zl*zl  +  OMEGAL)); 

tauh=0.0691*(1.0-Y)*0MEGAB*H; 

tau= (2 . 0/3 . 0) * (tauh/OMEGAM) * 

(sqrt(1.0-0MEGAM+0MEGAM*zl*zl*zl)-1.0); 

expt=exp(-2.0*tau) ; 

lamaf ac=(8 . 25e-7) *zl*zl ; 

return  £Lnswer*dxtodz*expt*zl*zl*lamaf ac/(extemp*extemp) ; 

} 

main  0 

{ 

int  t; 

double  ce,  ceec,  cee2,  cu,  ccl,  en,  ccon, 

constants,  11,  zl,  zzl,  z; 
double  qromo(double  (*func) (double) ,  double  a,  double  b, 

double  (*choose) (double (*) (double) ,  double,  double,  int)); 
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cu=H*100 . 0/ (3 . 0856el9*3 . 0856el9) ; 

/*  dimensional  and  unit  conversion  constants  */ 
ccl=XE*XE*OMEGAB*OMEGAB  *9 . 0*H*H* 10000 . 0*ST*ST 
/(64.0*PI*PI*G*G*MP*MP);   /*  electron 
density  squeired  times  cross  section  squcired 
/Hnaught  squeired.  */ 

cn=4.0*PI*PI;  /*  correction  for  difference  in 
normalization  between  Fry  and  me  */ 

ccon=H*H*10000.0/(8.0*PI*PI) ; 
constants=cu*ccl*cn*ccon/ (SOL+SOL) ; 

for  (t=0;t<=TMAX;t++)  { 

11=1 .0+(t*LLMAX/TMAX) ; 

l=pow(10.0,  11); 

ce=qromo(z integrand,  1.0,  21.0,  midpnt) ; 

ceec=constants*ce ; 

cee2=ceec*l* (1+1 . 0) / (2 . 0*PI) ; 

printf("y,d  '/.eXn",  1,  cee2) ; 


'^ 
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